Functional Data Analysis of London Bike Sharing Data
WJ
2024-12-02
Initial Setups
Introduction
This analysis includes data preparation, exploratory data analysis, Fourier smoothing, FPCA, FANOVA, functional regression, depth analysis, and clustering. All plot sets will be displayed in grids for easy comparison.
Data Preparation and Initial Checking
bike_data <- read.csv("london_merged.csv")
bike_data$timestamp <- as.POSIXct(bike_data$timestamp, format = "%Y-%m-%d %H:%M:%S")
str(bike_data)## 'data.frame': 17414 obs. of 10 variables:
## $ timestamp : POSIXct, format: "2015-01-04 00:00:00" "2015-01-04 01:00:00" ...
## $ cnt : int 182 138 134 72 47 46 51 75 131 301 ...
## $ t1 : num 3 3 2.5 2 2 2 1 1 1.5 2 ...
## $ t2 : num 2 2.5 2.5 2 0 2 -1 -1 -1 -0.5 ...
## $ hum : num 93 93 96.5 100 93 93 100 100 96.5 100 ...
## $ wind_speed : num 6 5 0 0 6.5 4 7 7 8 9 ...
## $ weather_code: num 3 1 1 1 1 1 4 4 4 3 ...
## $ is_holiday : num 0 0 0 0 0 0 0 0 0 0 ...
## $ is_weekend : num 1 1 1 1 1 1 1 1 1 1 ...
## $ season : num 3 3 3 3 3 3 3 3 3 3 ...## timestamp cnt t1 t2 hum wind_speed weather_code is_holiday
## 1 2015-01-04 00:00:00 182 3.0 2.0 93.0 6.0 3 0
## 2 2015-01-04 01:00:00 138 3.0 2.5 93.0 5.0 1 0
## 3 2015-01-04 02:00:00 134 2.5 2.5 96.5 0.0 1 0
## 4 2015-01-04 03:00:00 72 2.0 2.0 100.0 0.0 1 0
## 5 2015-01-04 04:00:00 47 2.0 0.0 93.0 6.5 1 0
## 6 2015-01-04 05:00:00 46 2.0 2.0 93.0 4.0 1 0
## is_weekend season
## 1 1 3
## 2 1 3
## 3 1 3
## 4 1 3
## 5 1 3
## 6 1 3##
## Dimensions: 17414 10##
## Missing values:## timestamp cnt t1 t2 hum wind_speed
## 0 0 0 0 0 0
## weather_code is_holiday is_weekend season
## 0 0 0 0##
## Summary:## timestamp cnt t1
## Min. :2015-01-04 00:00:00.00 Min. : 0 Min. :-1.50
## 1st Qu.:2015-07-04 20:15:00.00 1st Qu.: 257 1st Qu.: 8.00
## Median :2016-01-03 15:30:00.00 Median : 844 Median :12.50
## Mean :2016-01-03 21:56:01.84 Mean :1143 Mean :12.47
## 3rd Qu.:2016-07-04 15:45:00.00 3rd Qu.:1672 3rd Qu.:16.00
## Max. :2017-01-03 23:00:00.00 Max. :7860 Max. :34.00
## t2 hum wind_speed weather_code
## Min. :-6.00 Min. : 20.50 Min. : 0.00 Min. : 1.000
## 1st Qu.: 6.00 1st Qu.: 63.00 1st Qu.:10.00 1st Qu.: 1.000
## Median :12.50 Median : 74.50 Median :15.00 Median : 2.000
## Mean :11.52 Mean : 72.32 Mean :15.91 Mean : 2.723
## 3rd Qu.:16.00 3rd Qu.: 83.00 3rd Qu.:20.50 3rd Qu.: 3.000
## Max. :34.00 Max. :100.00 Max. :56.50 Max. :26.000
## is_holiday is_weekend season
## Min. :0.00000 Min. :0.0000 Min. :0.000
## 1st Qu.:0.00000 1st Qu.:0.0000 1st Qu.:0.000
## Median :0.00000 Median :0.0000 Median :1.000
## Mean :0.02205 Mean :0.2854 Mean :1.492
## 3rd Qu.:0.00000 3rd Qu.:1.0000 3rd Qu.:2.000
## Max. :1.00000 Max. :1.0000 Max. :3.000Feature Engineering
bike_data <- bike_data %>%
mutate(
hour = hour(timestamp),
day = as.Date(timestamp),
month = month(timestamp, label = TRUE),
week = week(timestamp),
year = year(timestamp),
day_of_week = wday(timestamp, label = TRUE),
is_weekend = ifelse(day_of_week %in% c("Sat", "Sun"), 1, 0),
day_of_month = mday(timestamp),
season_label = factor(season, levels = 0:3, labels = c("Spring", "Summer", "Fall", "Winter")),
weather_label = factor(
weather_code,
levels = c(1, 2, 3, 4, 7, 10, 26, 94),
labels = c("Clear", "Scattered Clouds", "Broken Clouds", "Cloudy",
"Rain", "Rain with Thunderstorm", "Snowfall", "Freezing Fog")
)
)
str(bike_data)## 'data.frame': 17414 obs. of 19 variables:
## $ timestamp : POSIXct, format: "2015-01-04 00:00:00" "2015-01-04 01:00:00" ...
## $ cnt : int 182 138 134 72 47 46 51 75 131 301 ...
## $ t1 : num 3 3 2.5 2 2 2 1 1 1.5 2 ...
## $ t2 : num 2 2.5 2.5 2 0 2 -1 -1 -1 -0.5 ...
## $ hum : num 93 93 96.5 100 93 93 100 100 96.5 100 ...
## $ wind_speed : num 6 5 0 0 6.5 4 7 7 8 9 ...
## $ weather_code : num 3 1 1 1 1 1 4 4 4 3 ...
## $ is_holiday : num 0 0 0 0 0 0 0 0 0 0 ...
## $ is_weekend : num 0 0 0 0 0 0 0 0 0 0 ...
## $ season : num 3 3 3 3 3 3 3 3 3 3 ...
## $ hour : int 0 1 2 3 4 5 6 7 8 9 ...
## $ day : Date, format: "2015-01-03" "2015-01-04" ...
## $ month : Ord.factor w/ 12 levels "1月"<"2月"<"3月"<..: 1 1 1 1 1 1 1 1 1 1 ...
## $ week : num 1 1 1 1 1 1 1 1 1 1 ...
## $ year : num 2015 2015 2015 2015 2015 ...
## $ day_of_week : Ord.factor w/ 7 levels "周日"<"周一"<..: 1 1 1 1 1 1 1 1 1 1 ...
## $ day_of_month : int 4 4 4 4 4 4 4 4 4 4 ...
## $ season_label : Factor w/ 4 levels "Spring","Summer",..: 4 4 4 4 4 4 4 4 4 4 ...
## $ weather_label: Factor w/ 8 levels "Clear","Scattered Clouds",..: 3 1 1 1 1 1 4 4 4 3 ...Create Daily Counts Matrix
daily_counts <- bike_data %>%
group_by(day, hour) %>%
summarise(cnt = mean(cnt, na.rm = TRUE), .groups = 'drop') %>%
pivot_wider(names_from = hour, values_from = cnt) %>%
ungroup() %>%
select(-day) %>%
as.matrix()
daily_counts[is.na(daily_counts)] <- 0
cat("\nDimensions of daily counts matrix:", dim(daily_counts), "\n")##
## Dimensions of daily counts matrix: 731 24Exploratory Plots
Time Series Patterns at Different Scales
daily_hourly <- bike_data %>%
group_by(day = as.Date(timestamp), hour) %>%
summarise(total_rentals = sum(cnt, na.rm = TRUE), .groups = "drop")
p1 <- ggplot(daily_hourly, aes(x = hour, y = total_rentals, group = day, color = as.factor(day))) +
geom_line(alpha = 0.6) + theme_minimal() + theme(legend.position = "none") +
labs(title = "Daily Rentals by Hour of Day", x = "Hour of Day", y = "Total Rentals")
weekly_daily <- bike_data %>%
group_by(week, day_of_week) %>%
summarise(total_rentals = sum(cnt, na.rm = TRUE), .groups = "drop")
p2 <- ggplot(weekly_daily, aes(x = day_of_week, y = total_rentals, group = week, color = as.factor(week))) +
geom_line(alpha = 0.6) + theme_minimal() + theme(legend.position = "none") +
labs(title = "Weekly Rentals by Day of Week", x = "Day of Week", y = "Total Rentals")
monthly_daily <- bike_data %>%
group_by(month, day_of_month) %>%
summarise(total_rentals = sum(cnt, na.rm = TRUE), .groups = "drop")
p3 <- ggplot(monthly_daily, aes(x = day_of_month, y = total_rentals, group = month, color = as.factor(month))) +
geom_line(alpha = 0.6) + theme_minimal() + theme(legend.position = "none") +
labs(title = "Monthly Rentals by Day of Month", x = "Day of Month", y = "Total Rentals")
yearly_monthly <- bike_data %>%
group_by(year, month) %>%
summarise(total_rentals = sum(cnt, na.rm = TRUE), .groups = "drop")
p4 <- ggplot(yearly_monthly, aes(x = month, y = total_rentals, group = year, color = as.factor(year))) +
geom_line(size = 1) + geom_point(size = 2) + theme_minimal() +
labs(title = "Yearly Rentals by Month of Year", x = "Month of Year", y = "Total Rentals")
# Arrange in a grid
grid.arrange(p1, p2, p3, p4, nrow=2, ncol=2)
Heatmaps
daily_hourly_heatmap <- bike_data %>%
group_by(day = as.Date(timestamp), hour) %>%
summarise(total_rentals = sum(cnt, na.rm = TRUE), .groups = 'drop')
p5 <- ggplot(daily_hourly_heatmap, aes(x = hour, y = day, fill = total_rentals)) +
geom_tile(color = "white") + scale_fill_gradient(low = "white", high = "blue") + theme_minimal() +
labs(title = "Daily-Hourly Bike Rentals", x = "Hour of Day", y = "Day")
weekly_daily_heatmap <- bike_data %>%
group_by(week, day_of_week) %>%
summarise(total_rentals = sum(cnt, na.rm = TRUE), .groups = 'drop')
p6 <- ggplot(weekly_daily_heatmap, aes(x = day_of_week, y = week, fill = total_rentals)) +
geom_tile(color = "white") + scale_fill_gradient(low = "white", high = "blue") + theme_minimal() +
labs(title = "Weekly-Daily Bike Rentals", x = "Day of Week", y = "Week")
monthly_daily_heatmap <- bike_data %>%
group_by(month, day_of_month) %>%
summarise(total_rentals = sum(cnt, na.rm = TRUE), .groups = 'drop')
p7 <- ggplot(monthly_daily_heatmap, aes(x = day_of_month, y = month, fill = total_rentals)) +
geom_tile(color = "white") + scale_fill_gradient(low = "white", high = "blue") + theme_minimal() +
labs(title = "Monthly-Daily Bike Rentals", x = "Day of Month", y = "Month")
yearly_monthly_heatmap <- bike_data %>%
group_by(year, month) %>%
summarise(total_rentals = sum(cnt, na.rm = TRUE), .groups = 'drop')
p8 <- ggplot(yearly_monthly_heatmap, aes(x = month, y = year, fill = total_rentals)) +
geom_tile(color = "white") + scale_fill_gradient(low = "white", high = "blue") + theme_minimal() +
labs(title = "Yearly-Monthly Bike Rentals", x = "Month", y = "Year")
grid.arrange(p5, p6, p7, p8, nrow=2, ncol=2)
Focused Heatmaps
weekly_hourly_focus <- bike_data %>%
group_by(day_of_week, hour) %>%
summarise(total_rentals = sum(cnt, na.rm = TRUE), .groups = 'drop')
p9 <- ggplot(weekly_hourly_focus, aes(x = hour, y = day_of_week, fill = total_rentals)) +
geom_tile(color = "white") + scale_fill_gradient(low = "white", high = "blue") + theme_minimal() +
labs(title = "Hourly Rentals by Day of Week", x = "Hour of Day", y = "Day of Week")
monthly_hourly_focus <- bike_data %>%
group_by(day_of_month, hour) %>%
summarise(total_rentals = sum(cnt, na.rm = TRUE), .groups = 'drop')
p10 <- ggplot(monthly_hourly_focus, aes(x = hour, y = day_of_month, fill = total_rentals)) +
geom_tile(color = "white") + scale_fill_gradient(low = "white", high = "blue") + theme_minimal() +
labs(title = "Hourly Rentals by Day of Month", x = "Hour of Day", y = "Day of Month")
yearly_hourly_focus <- bike_data %>%
group_by(month, hour) %>%
summarise(total_rentals = sum(cnt, na.rm = TRUE), .groups = 'drop')
p11 <- ggplot(yearly_hourly_focus, aes(x = hour, y = month, fill = total_rentals)) +
geom_tile(color = "white") + scale_fill_gradient(low = "white", high = "blue") + theme_minimal() +
labs(title = "Hourly Rentals by Month of Year", x = "Hour of Day", y = "Month")
grid.arrange(p9, p10, p11, ncol = 2)
Hourly Trends by Season, Month, and Holiday
season_hourly <- bike_data %>%
group_by(season_label, hour) %>%
summarise(mean_rentals = mean(cnt, na.rm = TRUE), .groups = 'drop')
p_season_grid <- ggplot(season_hourly, aes(x = hour, y = mean_rentals, color = season_label)) +
geom_line(size = 1) + facet_wrap(~season_label, ncol = 2) + theme_minimal() + theme(legend.position = "none") +
labs(title = "Hourly Trends by Season", x = "Hour of Day", y = "Mean Rentals")
month_hourly <- bike_data %>%
group_by(month, hour) %>%
summarise(mean_rentals = mean(cnt, na.rm = TRUE), .groups = 'drop')
p_month_grid <- ggplot(month_hourly, aes(x = hour, y = mean_rentals, color = month)) +
geom_line(size = 1) + facet_wrap(~month, ncol = 4) + theme_minimal() + theme(legend.position = "none") +
labs(title = "Hourly Trends by Month", x = "Hour of Day", y = "Mean Rentals")
holiday_hourly <- bike_data %>%
group_by(is_holiday, hour) %>%
summarise(mean_rentals = mean(cnt, na.rm = TRUE), .groups = 'drop')
p_holiday_grid <- ggplot(holiday_hourly, aes(x = hour, y = mean_rentals, color = as.factor(is_holiday))) +
geom_line(size = 1) + facet_wrap(~is_holiday, labeller = labeller(is_holiday = c("0" = "Non-Holiday", "1" = "Holiday"))) +
theme_minimal() + theme(legend.position = "none") +
labs(title = "Hourly Trends: Holiday vs Non-Holiday", x = "Hour of Day", y = "Mean Rentals")
# These are facetted plots themselves, but just show them one by one or arrange them together.
grid.arrange(p_season_grid, p_month_grid, p_holiday_grid, nrow = 2, ncol=2)
Functional Data Analysis
Fourier Smoothing with Dynamic Parameter Selection
hourly_data <- bike_data %>%
group_by(day, hour) %>%
summarise(
cnt = mean(cnt, na.rm = TRUE),
t1 = mean(t1, na.rm = TRUE),
t2 = mean(t2, na.rm = TRUE),
hum = mean(hum, na.rm = TRUE),
wind_speed = mean(wind_speed, na.rm = TRUE),
.groups = 'drop'
) %>%
pivot_wider(names_from = hour, values_from = c(cnt, t1, t2, hum, wind_speed), values_fill = 0)
variables <- c("cnt", "t1", "t2", "hum", "wind_speed")
basis <- create.fourier.basis(rangeval = c(0, 24), nbasis = nbasis)
results <- list()
for (var in variables) {
variable_data <- hourly_data %>%
select(starts_with(var)) %>%
as.matrix()
smoothed_variable <- matrix(0, nrow = nrow(variable_data), ncol = ncol(variable_data))
for (i in 1:nrow(variable_data)) {
y <- variable_data[i, ]
lambda_range <- 10^seq(-6, 0, by = 0.25)
gcv_values <- sapply(lambda_range, function(lam) {
fdParObj <- fdPar(basis, Lfdobj = 2, lambda = lam)
smooth_obj <- smooth.basis(seq(0, 24, length.out = length(y)), y, fdParObj)
smooth_obj$gcv
})
optimal_lambda <- lambda_range[which.min(gcv_values)]
fdParObj <- fdPar(basis, Lfdobj = 2, lambda = optimal_lambda)
smooth_obj <- smooth.basis(seq(0, 24, length.out = length(y)), y, fdParObj)
smoothed_variable[i, ] <- eval.fd(seq(0, 24, length.out = length(y)), smooth_obj$fd)
}
results[[var]] <- list(original = variable_data, smoothed = smoothed_variable)
}
# Prepare a list to store all plots
plot_list_os <- list()
for (var in variables) {
original_data <- results[[var]]$original
smoothed_data <- results[[var]]$smoothed
original_df <- as.data.frame(t(original_data))
smoothed_df <- as.data.frame(t(smoothed_data))
colnames(original_df) <- paste0("Day_", 1:nrow(original_data))
colnames(smoothed_df) <- paste0("Day_", 1:nrow(smoothed_data))
original_df$Hour <- seq(0, 23)
smoothed_df$Hour <- seq(0, 23)
original_long <- pivot_longer(original_df, -Hour, names_to = "Day", values_to = "Value")
smoothed_long <- pivot_longer(smoothed_df, -Hour, names_to = "Day", values_to = "Value")
p_o <- ggplot(original_long, aes(x = Hour, y = Value, color = Day)) +
geom_line(alpha = 0.5) + theme_minimal() + theme(legend.position = "none") +
labs(title = paste("Original", var), x = "Hour of Day", y = var)
p_s <- ggplot(smoothed_long, aes(x = Hour, y = Value, color = Day)) +
geom_line(alpha = 0.5) + theme_minimal() + theme(legend.position = "none") +
labs(title = paste("Smoothed", var), x = "Hour of Day", y = var)
# Add to the list
plot_list_os[[paste0(var, "_original")]] <- p_o
plot_list_os[[paste0(var, "_smoothed")]] <- p_s
}
# Arrange all 10 plots (2 for each of the 5 variables) in a 2x5 grid
grid.arrange(grobs = plot_list_os, nrow = 3, ncol = 4)
Functional Principal Component Analysis (FPCA)
pca_results <- list()
fpca_plots_var <- list()
fpca_plots_amp <- list() # store amplitude plots separately
for (var in variables) {
smoothed_data <- results[[var]]$smoothed
smoothed_fdObj <- Data2fd(
argvals = seq(0, 24, length.out = ncol(smoothed_data)),
y = t(smoothed_data),
basisobj = basis
)
fPCA <- pca.fd(smoothed_fdObj, nharm = max_nharm)
total_variance_threshold <- 0.90
var_prop <- fPCA$varprop
cum_var <- cumsum(var_prop)
dynamic_nharm <- which(cum_var >= total_variance_threshold)[1]
cat("Variable:", var, "\n")
cat("Number of PCs chosen dynamically:", dynamic_nharm, "\n")
cat("Cumulative variance explained:", cum_var[dynamic_nharm], "\n\n")
pca_results[[var]] <- list(
fPCA = fPCA,
dynamic_nharm = dynamic_nharm,
variance_explained = cum_var[dynamic_nharm]
)
# Create variance explained plot with cumulative line
variance_df <- data.frame(PC = 1:max_nharm, Variance = var_prop, Cumulative = cum_var)
p_var <- ggplot(variance_df, aes(x = PC)) +
geom_bar(aes(y = Variance), stat = "identity", fill = "steelblue", alpha = 0.7) +
geom_line(aes(y = Cumulative), color = "red", size = 1, group = 1) +
geom_point(aes(y = Cumulative), color = "red", size = 2) +
labs(title = paste("Variance Explained for", var),
x = "Principal Component",
y = "Proportion of Variance") +
theme_minimal() +
# Annotate chosen number of PCs
annotate("text", x = dynamic_nharm, y = cum_var[dynamic_nharm] + 0.05,
label = paste0("Chosen PCs: ", dynamic_nharm, "\nCumulative: ", round(cum_var[dynamic_nharm]*100, 2), "%"),
color = "black", size = 4, hjust = 0.5)
fpca_plots_var[[var]] <- p_var
# Varimax rotation for later plotting
fPCA_varimax <- varmx.pca.fd(fPCA)
# Store these objects for amplitude plots after loop
pca_results[[var]]$fPCA_varimax <- fPCA_varimax
}## Variable: cnt
## Number of PCs chosen dynamically: 2
## Cumulative variance explained: 0.9323043
##
## Variable: t1
## Number of PCs chosen dynamically: 1
## Cumulative variance explained: 0.91253
##
## Variable: t2
## Number of PCs chosen dynamically: 1
## Cumulative variance explained: 0.9200775
##
## Variable: hum
## Number of PCs chosen dynamically: 4
## Cumulative variance explained: 0.9007801
##
## Variable: wind_speed
## Number of PCs chosen dynamically: 4
## Cumulative variance explained: 0.9184972# Print variance explained plots in a grid
# Adjust rows/cols as needed. For example, if 5 variables, 2 rows by 3 cols:
all_plots_count <- length(fpca_plots_var)
dummy_plot <- ggplot() + theme_void() + labs(title=" ")
if (all_plots_count < 6) {
# Add dummy plots if fewer than 6
for (i in seq(all_plots_count+1, 6)) {
fpca_plots_var[[paste0("dummy", i)]] <- dummy_plot
}
}
grid.arrange(grobs = fpca_plots_var, nrow = 2, ncol = 3)
# For harmonic plots (amplitude):
# use matplot and legend for clear identification of PCs.
par(mfrow=c(2,3))
for (var in names(pca_results)) {
fPCA_obj <- pca_results[[var]]$fPCA
dynamic_nharm <- pca_results[[var]]$dynamic_nharm
time_grid <- seq(0, 24, length.out = 100)
# Evaluate harmonics for the selected PCs
harm_matrix <- eval.fd(time_grid, fPCA_obj$harmonics[1:dynamic_nharm])
matplot(time_grid, harm_matrix, type='l', lty=1, col=rainbow(dynamic_nharm),
xlab="Hour of Day", ylab="Amplitude",
main=paste("Principal Component Functions for", var))
legend("topright", legend=paste("PC", 1:dynamic_nharm),
col=rainbow(dynamic_nharm), lty=1, cex=0.8)
}
# If you have fewer than 6 variables, fill the last cell with a blank plot if needed
plot.new()
par(mfrow=c(2,3))
for (var in names(pca_results)) {
fPCA_obj <- pca_results[[var]]$fPCA
dynamic_nharm <- pca_results[[var]]$dynamic_nharm
fPCA_varimax <- pca_results[[var]]$fPCA_varimax
time_grid <- seq(0, 24, length.out = 100)
harm_matrix_varimax <- eval.fd(time_grid, fPCA_varimax$harmonics[1:dynamic_nharm])
matplot(time_grid, harm_matrix_varimax, type='l', lty=1, col=rainbow(dynamic_nharm),
xlab="Hour of Day", ylab="Amplitude",
main=paste("PCs After VARIMAX for", var))
legend("topright", legend=paste("PC", 1:dynamic_nharm),
col=rainbow(dynamic_nharm), lty=1, cex=0.8)
}
plot.new()
Functional ANOVA (FANOVA)
hourly_data_cnt <- bike_data %>%
group_by(day, hour) %>%
summarise(cnt = mean(cnt, na.rm = TRUE), .groups = 'drop') %>%
pivot_wider(names_from = hour, values_from = cnt, values_fill = 0)
basis_fanova <- create.fourier.basis(rangeval = c(0, 24), nbasis = nbasis)
smoothed_data_cnt <- t(apply(hourly_data_cnt[-1], 1, function(y) {
smooth_basis <- smooth.basis(seq(0, 23, length.out = length(y)), y, fdPar(basis_fanova))
eval.fd(seq(0, 23, length.out = length(y)), smooth_basis$fd)
}))
smoothed_data_t_cnt <- t(smoothed_data_cnt)
cnt_fd <- Data2fd(argvals = seq(0, 23, length.out = nrow(smoothed_data_t_cnt)),
y = smoothed_data_t_cnt, basisobj = basis_fanova)
unique_days <- unique(bike_data$day)
season_labels <- bike_data %>%
group_by(day) %>%
summarise(season = first(season_label)) %>%
arrange(day) %>%
pull(season)
holiday_labels <- bike_data %>%
group_by(day) %>%
summarise(is_holiday = first(is_holiday)) %>%
arrange(day) %>%
pull(is_holiday)
weather_labels <- bike_data %>%
group_by(day) %>%
summarise(weather = first(weather_label)) %>%
arrange(day) %>%
pull(weather)
run_fanova <- function(fd_obj, group_factor, group_name, basis_obj) {
groups <- factor(group_factor)
zmat <- model.matrix(~ groups - 1)
n_beta_basis <- nbasis
beta_basis_fanova <- create.fourier.basis(rangeval = c(0, 24), nbasis = n_beta_basis)
harm_accel_Lfd <- vec2Lfd(c((2 * pi / 24)^2, 0, 1), rangeval = c(0, 24))
beta_fd <- fd(matrix(0, n_beta_basis, ncol(zmat)), beta_basis_fanova)
beta_fdPar <- fdPar(beta_fd, harm_accel_Lfd, lambda = 1e-5)
beta_list_fanova <- lapply(1:ncol(zmat), function(i) beta_fdPar)
xfdlist_fanova <- lapply(1:ncol(zmat), function(i) zmat[, i])
fRegress_res <- fRegress(fd_obj, xfdlist_fanova, beta_list_fanova)
predicted_values_fd <- fRegress_res$yhatfdobj
observed_values_fd <- fd_obj
residuals_fd <- observed_values_fd$coefs - predicted_values_fd$coefs
SSE <- sum(residuals_fd^2)
SST <- sum((as.vector(observed_values_fd$coefs) - mean(as.vector(observed_values_fd$coefs)))^2)
R_squared <- 1 - SSE / SST
list(Group = group_name, SSE = SSE, SST = SST, R_squared = R_squared, BetaEst = fRegress_res$betaestlist)
}
fanova_season <- run_fanova(cnt_fd, season_labels, "Season", basis_fanova)
fanova_holiday <- run_fanova(cnt_fd, holiday_labels, "Holiday", basis_fanova)
fanova_weather <- run_fanova(cnt_fd, weather_labels, "Weather", basis_fanova)
cat("\nFANOVA Results:\n")##
## FANOVA Results:fanova_results_table <- data.frame(
Grouping = c(fanova_season$Group, fanova_holiday$Group, fanova_weather$Group),
SSE = c(fanova_season$SSE, fanova_holiday$SSE, fanova_weather$SSE),
SST = c(fanova_season$SST, fanova_holiday$SST, fanova_weather$SST),
R_squared = c(fanova_season$R_squared, fanova_holiday$R_squared, fanova_weather$R_squared)
)
print(fanova_results_table)## Grouping SSE SST R_squared
## 1 Season 6812857878 43223274248 0.8423799
## 2 Holiday 8037902084 43223274248 0.8140376
## 3 Weather 8027805070 43223274248 0.8142712# Plot Beta Functions in a grid for each factor set
# use par(mfrow) to arrange them
beta_plots_list <- list()
par(mfrow=c(2,3))
for (i in seq_along(fanova_season$BetaEst)) {
plot(fanova_season$BetaEst[[i]]$fd, main = paste("Season:", colnames(model.matrix(~season_labels - 1))[i]),
xlab = "Hour of Day", ylab = "Effect")
}
par(mfrow=c(2,3))
for (i in seq_along(fanova_holiday$BetaEst)) {
plot(fanova_holiday$BetaEst[[i]]$fd, main = paste("Holiday:", colnames(model.matrix(~holiday_labels - 1))[i]),
xlab = "Hour of Day", ylab = "Effect")
}
par(mfrow=c(2,3))
for (i in seq_along(fanova_weather$BetaEst)) {
plot(fanova_weather$BetaEst[[i]]$fd, main = paste("Weather:", colnames(model.matrix(~weather_labels - 1))[i]),
xlab = "Hour of Day", ylab = "Effect")
}
Functional Regression
predictor_data <- bike_data %>%
group_by(day, hour) %>%
summarise(
cnt = mean(cnt, na.rm = TRUE),
t2 = mean(t2, na.rm = TRUE),
hum = mean(hum, na.rm = TRUE),
wind_speed = mean(wind_speed, na.rm = TRUE),
holiday = first(is_holiday),
.groups = 'drop'
) %>%
arrange(day, hour)
holiday_daily <- predictor_data %>% distinct(day, holiday) %>% arrange(day)
predictor_wide <- predictor_data %>%
pivot_wider(names_from = hour,
values_from = c(cnt, t2, hum, wind_speed),
values_fill = 0) %>%
arrange(day)
holiday_vec <- holiday_daily$holiday
cnt_mat <- predictor_wide %>% select(starts_with("cnt_")) %>% as.matrix()
t2_mat <- predictor_wide %>% select(starts_with("t2_")) %>% as.matrix()
hum_mat <- predictor_wide %>% select(starts_with("hum_")) %>% as.matrix()
wind_mat<- predictor_wide %>% select(starts_with("wind_speed_")) %>% as.matrix()
make_fd <- function(data_matrix, basis) {
smooth_data_t <- t(apply(data_matrix, 1, function(y) {
sb <- smooth.basis(seq(0, 23, length.out = length(y)), y, fdPar(basis))
eval.fd(seq(0, 23, length.out = length(y)), sb$fd)
}))
Data2fd(seq(0, 23, length.out = ncol(smooth_data_t)), t(smooth_data_t), basis)
}
basis_reg <- create.fourier.basis(rangeval = c(0, 24), nbasis = nbasis)
cnt_fd_reg <- make_fd(cnt_mat, basis_reg)
t2_fd <- make_fd(t2_mat, basis_reg)
hum_fd <- make_fd(hum_mat, basis_reg)
wind_fd <- make_fd(wind_mat, basis_reg)
xfdlist_with_const <- list(
const = rep(1, nrow(cnt_mat)),
t2 = t2_fd,
hum = hum_fd,
wind = wind_fd,
holiday = holiday_vec
)
xfdlist_no_const <- list(
t2 = t2_fd,
hum = hum_fd,
wind = wind_fd,
holiday = holiday_vec
)
beta0_fdPar <- fdPar(create.constant.basis(c(0,24)))
beta_fdPar <- fdPar(basis_reg)
beta_list_with_const <- list(
beta0_fdPar,
beta_fdPar,
beta_fdPar,
beta_fdPar,
beta0_fdPar
)
beta_list_no_const <- list(
beta_fdPar,
beta_fdPar,
beta_fdPar,
beta0_fdPar
)
fRegress_with_const <- fRegress(cnt_fd_reg, xfdlist_with_const, beta_list_with_const)
pred_with_const <- fRegress_with_const$yhatfdobj
res_with_const <- cnt_fd_reg$coefs - pred_with_const$coefs
SSE_with_const <- sum(res_with_const^2)
SST_with_const <- sum((as.vector(cnt_fd_reg$coefs) - mean(as.vector(cnt_fd_reg$coefs)))^2)
R2_with_const <- 1 - SSE_with_const/SST_with_const
fRegress_no_const <- fRegress(cnt_fd_reg, xfdlist_no_const, beta_list_no_const)
pred_no_const <- fRegress_no_const$yhatfdobj
res_no_const <- cnt_fd_reg$coefs - pred_no_const$coefs
SSE_no_const <- sum(res_no_const^2)
R2_no_const <- 1 - SSE_no_const/SST_with_const
regression_results <- data.frame(
Model = c("With Constant", "Without Constant"),
SSE = c(SSE_with_const, SSE_no_const),
SST = c(SST_with_const, SST_with_const),
R_squared = c(R2_with_const, R2_no_const)
)
cat("\nFunctional Regression Metrics:\n")##
## Functional Regression Metrics:## Model SSE SST R_squared
## 1 With Constant 6634794378 43223892203 0.8465017
## 2 Without Constant 6999049455 43223892203 0.8380745Model Diagnostics and Additional Metrics
# Residual Analysis for Model with Constant
predicted_values_with_const <- pred_with_const
residuals_with_const <- cnt_fd_reg$coefs - predicted_values_with_const$coefs
# Residual Analysis for Model without Constant
predicted_values_no_const <- pred_no_const
residuals_no_const <- cnt_fd_reg$coefs - predicted_values_no_const$coefs
# Function to calculate metrics
calculate_metrics <- function(observed, predicted) {
observed_vals <- as.vector(observed)
predicted_vals <- as.vector(predicted)
MSE <- mean((observed_vals - predicted_vals)^2)
RMSE <- sqrt(MSE)
MAE <- mean(abs(observed_vals - predicted_vals))
MAPE <- mean(abs((observed_vals - predicted_vals) / (observed_vals + 1e-6))) * 100
list(MSE = MSE, RMSE = RMSE, MAE = MAE, MAPE = MAPE)
}
# Calculate metrics for both models
metrics_with_const <- calculate_metrics(cnt_fd_reg$coefs, predicted_values_with_const$coefs)
metrics_no_const <- calculate_metrics(cnt_fd_reg$coefs, predicted_values_no_const$coefs)
cat("\nAdditional Model Evaluation Metrics (With Constant):\n")##
## Additional Model Evaluation Metrics (With Constant):## MSE: 379066.1## RMSE: 615.6835## MAE: 401.1913## MAPE: 11677569 %##
## Additional Model Evaluation Metrics (Without Constant):## MSE: 399877.1## RMSE: 632.3584## MAE: 405.4448## MAPE: 11840360 %# Plot Predicted vs Observed for each model
num_days <- ncol(cnt_fd_reg$coefs)
rows <- 4
cols <- 6
par(mfrow = c(rows, cols), mar = c(2, 2, 2, 1))
# With Constant
for (i in 1:num_days) {
plot(predicted_values_with_const$coefs[, i], cnt_fd_reg$coefs[, i],
main = paste("With Const: Day", i), xlab = "Predicted", ylab = "Observed",
pch = 16, cex = 0.5)
abline(a = 0, b = 1, col = "red")
}
for (i in 1:num_days) {
plot(predicted_values_no_const$coefs[, i], cnt_fd_reg$coefs[, i],
main = paste("No Const: Day", i), xlab = "Predicted", ylab = "Observed",
pch = 16, cex = 0.5)
abline(a = 0, b = 1, col = "red")
}
# Residual Plot Comparison
par(mfrow = c(1, 2))
plot(as.vector(residuals_with_const), pch = 16,
main = "Residuals Plot (With Const)", xlab = "Observation Index", ylab = "Residuals")
abline(h = 0, col = "blue")
plot(as.vector(residuals_no_const), pch = 16,
main = "Residuals Plot (No Const)", xlab = "Observation Index", ylab = "Residuals")
abline(h = 0, col = "blue")
par(mfrow = c(1, 1))
# Residual Histogram Comparison
library(ggplot2)
library(gridExtra)
# With Constant
h_res_with_const <- ggplot(data.frame(residuals = as.vector(residuals_with_const)), aes(x = residuals)) +
geom_histogram(binwidth = diff(range(as.vector(residuals_with_const)))/30, fill = "skyblue", color = "white") +
labs(title = "Histogram of Residuals (With Const)", x = "Residual", y = "Frequency") +
theme_minimal()
# Without Constant
h_res_no_const <- ggplot(data.frame(residuals = as.vector(residuals_no_const)), aes(x = residuals)) +
geom_histogram(binwidth = diff(range(as.vector(residuals_no_const)))/30, fill = "orange", color = "white") +
labs(title = "Histogram of Residuals (No Const)", x = "Residual", y = "Frequency") +
theme_minimal()
# Combine Histograms
grid.arrange(h_res_with_const, h_res_no_const, ncol = 2)
# Q-Q Plot Comparison
qq_res_with_const <- ggplot(data.frame(residuals = as.vector(residuals_with_const)), aes(sample = residuals)) +
stat_qq() + stat_qq_line(col = "red") +
labs(title = "Q-Q Plot of Residuals (With Const)", x = "Theoretical Quantiles", y = "Sample Quantiles") +
theme_minimal()
qq_res_no_const <- ggplot(data.frame(residuals = as.vector(residuals_no_const)), aes(sample = residuals)) +
stat_qq() + stat_qq_line(col = "red") +
labs(title = "Q-Q Plot of Residuals (No Const)", x = "Theoretical Quantiles", y = "Sample Quantiles") +
theme_minimal()
# Combine Q-Q Plots
grid.arrange(qq_res_with_const, qq_res_no_const, ncol = 2)
# Print coefficients for model with constant
cat("Functional Coefficients for Model with Constant:\n")## Functional Coefficients for Model with Constant:for (i in seq_along(fRegress_with_const$betaestlist)) {
print(fRegress_with_const$betaestlist[[i]])
}## $fd
## $coefs
## [,1]
## [1,] 607.1884
##
## $basis
## $call
## basisfd(type = type, rangeval = rangeval, nbasis = nbasis, params = params,
## dropind = dropind, quadvals = quadvals, values = values,
## basisvalues = basisvalues)
##
## $type
## [1] "const"
##
## $rangeval
## [1] 0 24
##
## $nbasis
## [1] 1
##
## $params
## [1] 0
##
## $dropind
## numeric(0)
##
## $quadvals
## numeric(0)
##
## $values
## list()
##
## $basisvalues
## list()
##
## $names
## [1] "const"
##
## attr(,"class")
## [1] "basisfd"
##
## $fdnames
## $fdnames[[1]]
## [1] "const"
##
## $fdnames[[2]]
## [1] "reps 1"
##
## $fdnames[[3]]
## [1] "values"
##
##
## attr(,"class")
## [1] "fd"
##
## $Lfd
## $call
## Lfd(nderiv = m, bwtlist = bwtlist)
##
## $nderiv
## [1] 0
##
## $bwtlist
## list()
##
## attr(,"class")
## [1] "Lfd"
##
## $lambda
## [1] 0
##
## $estimate
## [1] TRUE
##
## $penmat
## NULL
##
## attr(,"class")
## [1] "fdPar"
## $fd
## $coefs
## [,1]
## [1,] 219.1297888
## [2,] -123.0421694
## [3,] -87.1106636
## [4,] -5.6456375
## [5,] -56.0833500
## [6,] 30.4388481
## [7,] 28.1709815
## [8,] 10.6214819
## [9,] 20.1578903
## [10,] -31.3152415
## [11,] -20.9235014
## [12,] 7.0108290
## [13,] 3.3606206
## [14,] 12.4618035
## [15,] 2.4138829
## [16,] -10.3874650
## [17,] -1.1396325
## [18,] -1.2341304
## [19,] 0.1400955
## [20,] 7.4402777
## [21,] -1.3474247
## [22,] -4.0202755
## [23,] -0.1976931
##
## $basis
## $call
## basisfd(type = type, rangeval = rangeval, nbasis = nbasis, params = params,
## dropind = dropind, quadvals = quadvals, values = values,
## basisvalues = basisvalues)
##
## $type
## [1] "fourier"
##
## $rangeval
## [1] 0 24
##
## $nbasis
## [1] 23
##
## $params
## [1] 24
##
## $dropind
## numeric(0)
##
## $quadvals
## NULL
##
## $values
## list()
##
## $basisvalues
## list()
##
## $names
## [1] "const" "sin1" "cos1" "sin2" "cos2" "sin3" "cos3" "sin4" "cos4"
## [10] "sin5" "cos5" "sin6" "cos6" "sin7" "cos7" "sin8" "cos8" "sin9"
## [19] "cos9" "sin10" "cos10" "sin11" "cos11"
##
## attr(,"class")
## [1] "basisfd"
##
## $fdnames
## $fdnames[[1]]
## [1] "const" "sin1" "cos1" "sin2" "cos2" "sin3" "cos3" "sin4" "cos4"
## [10] "sin5" "cos5" "sin6" "cos6" "sin7" "cos7" "sin8" "cos8" "sin9"
## [19] "cos9" "sin10" "cos10" "sin11" "cos11"
##
## $fdnames[[2]]
## [1] "reps 1"
##
## $fdnames[[3]]
## [1] "values"
##
##
## attr(,"class")
## [1] "fd"
##
## $Lfd
## $call
## Lfd(nderiv = m, bwtlist = bwtlist)
##
## $nderiv
## [1] 3
##
## $bwtlist
## $bwtlist[[1]]
## $coefs
## reps 1
## time 0
##
## $basis
## $call
## basisfd(type = type, rangeval = rangeval, nbasis = nbasis, params = params,
## dropind = dropind, quadvals = quadvals, values = values,
## basisvalues = basisvalues)
##
## $type
## [1] "const"
##
## $rangeval
## [1] 0 24
##
## $nbasis
## [1] 1
##
## $params
## [1] 0
##
## $dropind
## numeric(0)
##
## $quadvals
## numeric(0)
##
## $values
## list()
##
## $basisvalues
## list()
##
## $names
## [1] "const"
##
## attr(,"class")
## [1] "basisfd"
##
## $fdnames
## $fdnames$args
## [1] "time"
##
## $fdnames$reps
## [1] "reps 1"
##
## $fdnames$funs
## [1] "values"
##
##
## attr(,"class")
## [1] "fd"
##
## $bwtlist[[2]]
## $coefs
## reps 1
## time 0.06853892
##
## $basis
## $call
## basisfd(type = type, rangeval = rangeval, nbasis = nbasis, params = params,
## dropind = dropind, quadvals = quadvals, values = values,
## basisvalues = basisvalues)
##
## $type
## [1] "const"
##
## $rangeval
## [1] 0 24
##
## $nbasis
## [1] 1
##
## $params
## [1] 0
##
## $dropind
## numeric(0)
##
## $quadvals
## numeric(0)
##
## $values
## list()
##
## $basisvalues
## list()
##
## $names
## [1] "const"
##
## attr(,"class")
## [1] "basisfd"
##
## $fdnames
## $fdnames$args
## [1] "time"
##
## $fdnames$reps
## [1] "reps 1"
##
## $fdnames$funs
## [1] "values"
##
##
## attr(,"class")
## [1] "fd"
##
## $bwtlist[[3]]
## $coefs
## reps 1
## time 0
##
## $basis
## $call
## basisfd(type = type, rangeval = rangeval, nbasis = nbasis, params = params,
## dropind = dropind, quadvals = quadvals, values = values,
## basisvalues = basisvalues)
##
## $type
## [1] "const"
##
## $rangeval
## [1] 0 24
##
## $nbasis
## [1] 1
##
## $params
## [1] 0
##
## $dropind
## numeric(0)
##
## $quadvals
## numeric(0)
##
## $values
## list()
##
## $basisvalues
## list()
##
## $names
## [1] "const"
##
## attr(,"class")
## [1] "basisfd"
##
## $fdnames
## $fdnames$args
## [1] "time"
##
## $fdnames$reps
## [1] "reps 1"
##
## $fdnames$funs
## [1] "values"
##
##
## attr(,"class")
## [1] "fd"
##
##
## attr(,"class")
## [1] "Lfd"
##
## $lambda
## [1] 0
##
## $estimate
## [1] TRUE
##
## $penmat
## NULL
##
## attr(,"class")
## [1] "fdPar"
## $fd
## $coefs
## [,1]
## [1,] -3.6069719
## [2,] 4.5750006
## [3,] -13.7915674
## [4,] -11.2513823
## [5,] -6.8994582
## [6,] 2.5543304
## [7,] 10.9331243
## [8,] 5.7464765
## [9,] -1.1491406
## [10,] -6.5861566
## [11,] -5.6625457
## [12,] 0.5124968
## [13,] 4.0843311
## [14,] 3.8204789
## [15,] -0.8065256
## [16,] -2.9942881
## [17,] -1.5323447
## [18,] 0.1970428
## [19,] 1.7250624
## [20,] 2.0545186
## [21,] -0.6718584
## [22,] -1.7180737
## [23,] -0.9728825
##
## $basis
## $call
## basisfd(type = type, rangeval = rangeval, nbasis = nbasis, params = params,
## dropind = dropind, quadvals = quadvals, values = values,
## basisvalues = basisvalues)
##
## $type
## [1] "fourier"
##
## $rangeval
## [1] 0 24
##
## $nbasis
## [1] 23
##
## $params
## [1] 24
##
## $dropind
## numeric(0)
##
## $quadvals
## NULL
##
## $values
## list()
##
## $basisvalues
## list()
##
## $names
## [1] "const" "sin1" "cos1" "sin2" "cos2" "sin3" "cos3" "sin4" "cos4"
## [10] "sin5" "cos5" "sin6" "cos6" "sin7" "cos7" "sin8" "cos8" "sin9"
## [19] "cos9" "sin10" "cos10" "sin11" "cos11"
##
## attr(,"class")
## [1] "basisfd"
##
## $fdnames
## $fdnames[[1]]
## [1] "const" "sin1" "cos1" "sin2" "cos2" "sin3" "cos3" "sin4" "cos4"
## [10] "sin5" "cos5" "sin6" "cos6" "sin7" "cos7" "sin8" "cos8" "sin9"
## [19] "cos9" "sin10" "cos10" "sin11" "cos11"
##
## $fdnames[[2]]
## [1] "reps 1"
##
## $fdnames[[3]]
## [1] "values"
##
##
## attr(,"class")
## [1] "fd"
##
## $Lfd
## $call
## Lfd(nderiv = m, bwtlist = bwtlist)
##
## $nderiv
## [1] 3
##
## $bwtlist
## $bwtlist[[1]]
## $coefs
## reps 1
## time 0
##
## $basis
## $call
## basisfd(type = type, rangeval = rangeval, nbasis = nbasis, params = params,
## dropind = dropind, quadvals = quadvals, values = values,
## basisvalues = basisvalues)
##
## $type
## [1] "const"
##
## $rangeval
## [1] 0 24
##
## $nbasis
## [1] 1
##
## $params
## [1] 0
##
## $dropind
## numeric(0)
##
## $quadvals
## numeric(0)
##
## $values
## list()
##
## $basisvalues
## list()
##
## $names
## [1] "const"
##
## attr(,"class")
## [1] "basisfd"
##
## $fdnames
## $fdnames$args
## [1] "time"
##
## $fdnames$reps
## [1] "reps 1"
##
## $fdnames$funs
## [1] "values"
##
##
## attr(,"class")
## [1] "fd"
##
## $bwtlist[[2]]
## $coefs
## reps 1
## time 0.06853892
##
## $basis
## $call
## basisfd(type = type, rangeval = rangeval, nbasis = nbasis, params = params,
## dropind = dropind, quadvals = quadvals, values = values,
## basisvalues = basisvalues)
##
## $type
## [1] "const"
##
## $rangeval
## [1] 0 24
##
## $nbasis
## [1] 1
##
## $params
## [1] 0
##
## $dropind
## numeric(0)
##
## $quadvals
## numeric(0)
##
## $values
## list()
##
## $basisvalues
## list()
##
## $names
## [1] "const"
##
## attr(,"class")
## [1] "basisfd"
##
## $fdnames
## $fdnames$args
## [1] "time"
##
## $fdnames$reps
## [1] "reps 1"
##
## $fdnames$funs
## [1] "values"
##
##
## attr(,"class")
## [1] "fd"
##
## $bwtlist[[3]]
## $coefs
## reps 1
## time 0
##
## $basis
## $call
## basisfd(type = type, rangeval = rangeval, nbasis = nbasis, params = params,
## dropind = dropind, quadvals = quadvals, values = values,
## basisvalues = basisvalues)
##
## $type
## [1] "const"
##
## $rangeval
## [1] 0 24
##
## $nbasis
## [1] 1
##
## $params
## [1] 0
##
## $dropind
## numeric(0)
##
## $quadvals
## numeric(0)
##
## $values
## list()
##
## $basisvalues
## list()
##
## $names
## [1] "const"
##
## attr(,"class")
## [1] "basisfd"
##
## $fdnames
## $fdnames$args
## [1] "time"
##
## $fdnames$reps
## [1] "reps 1"
##
## $fdnames$funs
## [1] "values"
##
##
## attr(,"class")
## [1] "fd"
##
##
## attr(,"class")
## [1] "Lfd"
##
## $lambda
## [1] 0
##
## $estimate
## [1] TRUE
##
## $penmat
## NULL
##
## attr(,"class")
## [1] "fdPar"
## $fd
## $coefs
## [,1]
## [1,] 2.6779189
## [2,] -13.8702742
## [3,] -11.3293288
## [4,] -0.6871929
## [5,] -13.2261821
## [6,] 9.4983305
## [7,] 6.4103265
## [8,] 1.1815567
## [9,] 3.0214120
## [10,] -7.8871530
## [11,] -4.2845055
## [12,] 3.3374414
## [13,] 1.5484137
## [14,] 2.8945936
## [15,] -0.6789299
## [16,] -3.5875487
## [17,] 0.7671949
## [18,] -0.6020117
## [19,] -0.1865963
## [20,] 2.4419199
## [21,] -1.1690623
## [22,] -2.5472327
## [23,] 0.5104993
##
## $basis
## $call
## basisfd(type = type, rangeval = rangeval, nbasis = nbasis, params = params,
## dropind = dropind, quadvals = quadvals, values = values,
## basisvalues = basisvalues)
##
## $type
## [1] "fourier"
##
## $rangeval
## [1] 0 24
##
## $nbasis
## [1] 23
##
## $params
## [1] 24
##
## $dropind
## numeric(0)
##
## $quadvals
## NULL
##
## $values
## list()
##
## $basisvalues
## list()
##
## $names
## [1] "const" "sin1" "cos1" "sin2" "cos2" "sin3" "cos3" "sin4" "cos4"
## [10] "sin5" "cos5" "sin6" "cos6" "sin7" "cos7" "sin8" "cos8" "sin9"
## [19] "cos9" "sin10" "cos10" "sin11" "cos11"
##
## attr(,"class")
## [1] "basisfd"
##
## $fdnames
## $fdnames[[1]]
## [1] "const" "sin1" "cos1" "sin2" "cos2" "sin3" "cos3" "sin4" "cos4"
## [10] "sin5" "cos5" "sin6" "cos6" "sin7" "cos7" "sin8" "cos8" "sin9"
## [19] "cos9" "sin10" "cos10" "sin11" "cos11"
##
## $fdnames[[2]]
## [1] "reps 1"
##
## $fdnames[[3]]
## [1] "values"
##
##
## attr(,"class")
## [1] "fd"
##
## $Lfd
## $call
## Lfd(nderiv = m, bwtlist = bwtlist)
##
## $nderiv
## [1] 3
##
## $bwtlist
## $bwtlist[[1]]
## $coefs
## reps 1
## time 0
##
## $basis
## $call
## basisfd(type = type, rangeval = rangeval, nbasis = nbasis, params = params,
## dropind = dropind, quadvals = quadvals, values = values,
## basisvalues = basisvalues)
##
## $type
## [1] "const"
##
## $rangeval
## [1] 0 24
##
## $nbasis
## [1] 1
##
## $params
## [1] 0
##
## $dropind
## numeric(0)
##
## $quadvals
## numeric(0)
##
## $values
## list()
##
## $basisvalues
## list()
##
## $names
## [1] "const"
##
## attr(,"class")
## [1] "basisfd"
##
## $fdnames
## $fdnames$args
## [1] "time"
##
## $fdnames$reps
## [1] "reps 1"
##
## $fdnames$funs
## [1] "values"
##
##
## attr(,"class")
## [1] "fd"
##
## $bwtlist[[2]]
## $coefs
## reps 1
## time 0.06853892
##
## $basis
## $call
## basisfd(type = type, rangeval = rangeval, nbasis = nbasis, params = params,
## dropind = dropind, quadvals = quadvals, values = values,
## basisvalues = basisvalues)
##
## $type
## [1] "const"
##
## $rangeval
## [1] 0 24
##
## $nbasis
## [1] 1
##
## $params
## [1] 0
##
## $dropind
## numeric(0)
##
## $quadvals
## numeric(0)
##
## $values
## list()
##
## $basisvalues
## list()
##
## $names
## [1] "const"
##
## attr(,"class")
## [1] "basisfd"
##
## $fdnames
## $fdnames$args
## [1] "time"
##
## $fdnames$reps
## [1] "reps 1"
##
## $fdnames$funs
## [1] "values"
##
##
## attr(,"class")
## [1] "fd"
##
## $bwtlist[[3]]
## $coefs
## reps 1
## time 0
##
## $basis
## $call
## basisfd(type = type, rangeval = rangeval, nbasis = nbasis, params = params,
## dropind = dropind, quadvals = quadvals, values = values,
## basisvalues = basisvalues)
##
## $type
## [1] "const"
##
## $rangeval
## [1] 0 24
##
## $nbasis
## [1] 1
##
## $params
## [1] 0
##
## $dropind
## numeric(0)
##
## $quadvals
## numeric(0)
##
## $values
## list()
##
## $basisvalues
## list()
##
## $names
## [1] "const"
##
## attr(,"class")
## [1] "basisfd"
##
## $fdnames
## $fdnames$args
## [1] "time"
##
## $fdnames$reps
## [1] "reps 1"
##
## $fdnames$funs
## [1] "values"
##
##
## attr(,"class")
## [1] "fd"
##
##
## attr(,"class")
## [1] "Lfd"
##
## $lambda
## [1] 0
##
## $estimate
## [1] TRUE
##
## $penmat
## NULL
##
## attr(,"class")
## [1] "fdPar"
## $fd
## $coefs
## [,1]
## [1,] -435.4812
##
## $basis
## $call
## basisfd(type = type, rangeval = rangeval, nbasis = nbasis, params = params,
## dropind = dropind, quadvals = quadvals, values = values,
## basisvalues = basisvalues)
##
## $type
## [1] "const"
##
## $rangeval
## [1] 0 24
##
## $nbasis
## [1] 1
##
## $params
## [1] 0
##
## $dropind
## numeric(0)
##
## $quadvals
## numeric(0)
##
## $values
## list()
##
## $basisvalues
## list()
##
## $names
## [1] "const"
##
## attr(,"class")
## [1] "basisfd"
##
## $fdnames
## $fdnames[[1]]
## [1] "const"
##
## $fdnames[[2]]
## [1] "reps 1"
##
## $fdnames[[3]]
## [1] "values"
##
##
## attr(,"class")
## [1] "fd"
##
## $Lfd
## $call
## Lfd(nderiv = m, bwtlist = bwtlist)
##
## $nderiv
## [1] 0
##
## $bwtlist
## list()
##
## attr(,"class")
## [1] "Lfd"
##
## $lambda
## [1] 0
##
## $estimate
## [1] TRUE
##
## $penmat
## NULL
##
## attr(,"class")
## [1] "fdPar"# Plot coefficients for model with constant
par(mfrow = c(2, 3)) # Adjust layout if needed
for (i in seq_along(fRegress_with_const$betaestlist)) {
plot(fRegress_with_const$betaestlist[[i]]$fd,
main = paste("Coefficient Function for Predictor", i))
}
par(mfrow = c(1, 1))
# Print coefficients for model without constant
cat("Functional Coefficients for Model without Constant:\n")## Functional Coefficients for Model without Constant:## $fd
## $coefs
## [,1]
## [1,] 250.1731429
## [2,] -132.8714995
## [3,] -96.1502158
## [4,] -5.9883118
## [5,] -54.2997405
## [6,] 31.5218777
## [7,] 27.8178968
## [8,] 10.9195800
## [9,] 19.7990896
## [10,] -31.3876431
## [11,] -21.1571878
## [12,] 7.1200533
## [13,] 3.1830232
## [14,] 12.5013915
## [15,] 2.1433094
## [16,] -10.3757027
## [17,] -1.3149212
## [18,] -1.2340220
## [19,] -0.1075665
## [20,] 7.4413216
## [21,] -1.6522911
## [22,] -3.9877769
## [23,] -0.3410381
##
## $basis
## $call
## basisfd(type = type, rangeval = rangeval, nbasis = nbasis, params = params,
## dropind = dropind, quadvals = quadvals, values = values,
## basisvalues = basisvalues)
##
## $type
## [1] "fourier"
##
## $rangeval
## [1] 0 24
##
## $nbasis
## [1] 23
##
## $params
## [1] 24
##
## $dropind
## numeric(0)
##
## $quadvals
## NULL
##
## $values
## list()
##
## $basisvalues
## list()
##
## $names
## [1] "const" "sin1" "cos1" "sin2" "cos2" "sin3" "cos3" "sin4" "cos4"
## [10] "sin5" "cos5" "sin6" "cos6" "sin7" "cos7" "sin8" "cos8" "sin9"
## [19] "cos9" "sin10" "cos10" "sin11" "cos11"
##
## attr(,"class")
## [1] "basisfd"
##
## $fdnames
## $fdnames[[1]]
## [1] "const" "sin1" "cos1" "sin2" "cos2" "sin3" "cos3" "sin4" "cos4"
## [10] "sin5" "cos5" "sin6" "cos6" "sin7" "cos7" "sin8" "cos8" "sin9"
## [19] "cos9" "sin10" "cos10" "sin11" "cos11"
##
## $fdnames[[2]]
## [1] "reps 1"
##
## $fdnames[[3]]
## [1] "values"
##
##
## attr(,"class")
## [1] "fd"
##
## $Lfd
## $call
## Lfd(nderiv = m, bwtlist = bwtlist)
##
## $nderiv
## [1] 3
##
## $bwtlist
## $bwtlist[[1]]
## $coefs
## reps 1
## time 0
##
## $basis
## $call
## basisfd(type = type, rangeval = rangeval, nbasis = nbasis, params = params,
## dropind = dropind, quadvals = quadvals, values = values,
## basisvalues = basisvalues)
##
## $type
## [1] "const"
##
## $rangeval
## [1] 0 24
##
## $nbasis
## [1] 1
##
## $params
## [1] 0
##
## $dropind
## numeric(0)
##
## $quadvals
## numeric(0)
##
## $values
## list()
##
## $basisvalues
## list()
##
## $names
## [1] "const"
##
## attr(,"class")
## [1] "basisfd"
##
## $fdnames
## $fdnames$args
## [1] "time"
##
## $fdnames$reps
## [1] "reps 1"
##
## $fdnames$funs
## [1] "values"
##
##
## attr(,"class")
## [1] "fd"
##
## $bwtlist[[2]]
## $coefs
## reps 1
## time 0.06853892
##
## $basis
## $call
## basisfd(type = type, rangeval = rangeval, nbasis = nbasis, params = params,
## dropind = dropind, quadvals = quadvals, values = values,
## basisvalues = basisvalues)
##
## $type
## [1] "const"
##
## $rangeval
## [1] 0 24
##
## $nbasis
## [1] 1
##
## $params
## [1] 0
##
## $dropind
## numeric(0)
##
## $quadvals
## numeric(0)
##
## $values
## list()
##
## $basisvalues
## list()
##
## $names
## [1] "const"
##
## attr(,"class")
## [1] "basisfd"
##
## $fdnames
## $fdnames$args
## [1] "time"
##
## $fdnames$reps
## [1] "reps 1"
##
## $fdnames$funs
## [1] "values"
##
##
## attr(,"class")
## [1] "fd"
##
## $bwtlist[[3]]
## $coefs
## reps 1
## time 0
##
## $basis
## $call
## basisfd(type = type, rangeval = rangeval, nbasis = nbasis, params = params,
## dropind = dropind, quadvals = quadvals, values = values,
## basisvalues = basisvalues)
##
## $type
## [1] "const"
##
## $rangeval
## [1] 0 24
##
## $nbasis
## [1] 1
##
## $params
## [1] 0
##
## $dropind
## numeric(0)
##
## $quadvals
## numeric(0)
##
## $values
## list()
##
## $basisvalues
## list()
##
## $names
## [1] "const"
##
## attr(,"class")
## [1] "basisfd"
##
## $fdnames
## $fdnames$args
## [1] "time"
##
## $fdnames$reps
## [1] "reps 1"
##
## $fdnames$funs
## [1] "values"
##
##
## attr(,"class")
## [1] "fd"
##
##
## attr(,"class")
## [1] "Lfd"
##
## $lambda
## [1] 0
##
## $estimate
## [1] TRUE
##
## $penmat
## NULL
##
## attr(,"class")
## [1] "fdPar"
## $fd
## $coefs
## [,1]
## [1,] 26.2432648
## [2,] 5.6575817
## [3,] -12.4284824
## [4,] -11.5196035
## [5,] -6.9890216
## [6,] 2.4806669
## [7,] 11.0919702
## [8,] 5.7609288
## [9,] -1.2195528
## [10,] -6.5481244
## [11,] -5.6382041
## [12,] 0.5230600
## [13,] 4.1208623
## [14,] 3.7868401
## [15,] -0.7491247
## [16,] -2.9872422
## [17,] -1.4972348
## [18,] 0.1348025
## [19,] 1.7490978
## [20,] 2.0408031
## [21,] -0.6390626
## [22,] -1.8183761
## [23,] -0.9534914
##
## $basis
## $call
## basisfd(type = type, rangeval = rangeval, nbasis = nbasis, params = params,
## dropind = dropind, quadvals = quadvals, values = values,
## basisvalues = basisvalues)
##
## $type
## [1] "fourier"
##
## $rangeval
## [1] 0 24
##
## $nbasis
## [1] 23
##
## $params
## [1] 24
##
## $dropind
## numeric(0)
##
## $quadvals
## NULL
##
## $values
## list()
##
## $basisvalues
## list()
##
## $names
## [1] "const" "sin1" "cos1" "sin2" "cos2" "sin3" "cos3" "sin4" "cos4"
## [10] "sin5" "cos5" "sin6" "cos6" "sin7" "cos7" "sin8" "cos8" "sin9"
## [19] "cos9" "sin10" "cos10" "sin11" "cos11"
##
## attr(,"class")
## [1] "basisfd"
##
## $fdnames
## $fdnames[[1]]
## [1] "const" "sin1" "cos1" "sin2" "cos2" "sin3" "cos3" "sin4" "cos4"
## [10] "sin5" "cos5" "sin6" "cos6" "sin7" "cos7" "sin8" "cos8" "sin9"
## [19] "cos9" "sin10" "cos10" "sin11" "cos11"
##
## $fdnames[[2]]
## [1] "reps 1"
##
## $fdnames[[3]]
## [1] "values"
##
##
## attr(,"class")
## [1] "fd"
##
## $Lfd
## $call
## Lfd(nderiv = m, bwtlist = bwtlist)
##
## $nderiv
## [1] 3
##
## $bwtlist
## $bwtlist[[1]]
## $coefs
## reps 1
## time 0
##
## $basis
## $call
## basisfd(type = type, rangeval = rangeval, nbasis = nbasis, params = params,
## dropind = dropind, quadvals = quadvals, values = values,
## basisvalues = basisvalues)
##
## $type
## [1] "const"
##
## $rangeval
## [1] 0 24
##
## $nbasis
## [1] 1
##
## $params
## [1] 0
##
## $dropind
## numeric(0)
##
## $quadvals
## numeric(0)
##
## $values
## list()
##
## $basisvalues
## list()
##
## $names
## [1] "const"
##
## attr(,"class")
## [1] "basisfd"
##
## $fdnames
## $fdnames$args
## [1] "time"
##
## $fdnames$reps
## [1] "reps 1"
##
## $fdnames$funs
## [1] "values"
##
##
## attr(,"class")
## [1] "fd"
##
## $bwtlist[[2]]
## $coefs
## reps 1
## time 0.06853892
##
## $basis
## $call
## basisfd(type = type, rangeval = rangeval, nbasis = nbasis, params = params,
## dropind = dropind, quadvals = quadvals, values = values,
## basisvalues = basisvalues)
##
## $type
## [1] "const"
##
## $rangeval
## [1] 0 24
##
## $nbasis
## [1] 1
##
## $params
## [1] 0
##
## $dropind
## numeric(0)
##
## $quadvals
## numeric(0)
##
## $values
## list()
##
## $basisvalues
## list()
##
## $names
## [1] "const"
##
## attr(,"class")
## [1] "basisfd"
##
## $fdnames
## $fdnames$args
## [1] "time"
##
## $fdnames$reps
## [1] "reps 1"
##
## $fdnames$funs
## [1] "values"
##
##
## attr(,"class")
## [1] "fd"
##
## $bwtlist[[3]]
## $coefs
## reps 1
## time 0
##
## $basis
## $call
## basisfd(type = type, rangeval = rangeval, nbasis = nbasis, params = params,
## dropind = dropind, quadvals = quadvals, values = values,
## basisvalues = basisvalues)
##
## $type
## [1] "const"
##
## $rangeval
## [1] 0 24
##
## $nbasis
## [1] 1
##
## $params
## [1] 0
##
## $dropind
## numeric(0)
##
## $quadvals
## numeric(0)
##
## $values
## list()
##
## $basisvalues
## list()
##
## $names
## [1] "const"
##
## attr(,"class")
## [1] "basisfd"
##
## $fdnames
## $fdnames$args
## [1] "time"
##
## $fdnames$reps
## [1] "reps 1"
##
## $fdnames$funs
## [1] "values"
##
##
## attr(,"class")
## [1] "fd"
##
##
## attr(,"class")
## [1] "Lfd"
##
## $lambda
## [1] 0
##
## $estimate
## [1] TRUE
##
## $penmat
## NULL
##
## attr(,"class")
## [1] "fdPar"
## $fd
## $coefs
## [,1]
## [1,] 2.372383e+01
## [2,] -1.671277e+01
## [3,] -1.446060e+01
## [4,] -6.558164e-04
## [5,] -1.341927e+01
## [6,] 9.710769e+00
## [7,] 6.028092e+00
## [8,] 8.100872e-01
## [9,] 3.396451e+00
## [10,] -8.032724e+00
## [11,] -4.043424e+00
## [12,] 3.155279e+00
## [13,] 1.641191e+00
## [14,] 2.948564e+00
## [15,] -6.900351e-01
## [16,] -3.684383e+00
## [17,] 7.872497e-01
## [18,] -5.420047e-01
## [19,] -2.019536e-02
## [20,] 2.410169e+00
## [21,] -1.025462e+00
## [22,] -2.536341e+00
## [23,] 5.642935e-01
##
## $basis
## $call
## basisfd(type = type, rangeval = rangeval, nbasis = nbasis, params = params,
## dropind = dropind, quadvals = quadvals, values = values,
## basisvalues = basisvalues)
##
## $type
## [1] "fourier"
##
## $rangeval
## [1] 0 24
##
## $nbasis
## [1] 23
##
## $params
## [1] 24
##
## $dropind
## numeric(0)
##
## $quadvals
## NULL
##
## $values
## list()
##
## $basisvalues
## list()
##
## $names
## [1] "const" "sin1" "cos1" "sin2" "cos2" "sin3" "cos3" "sin4" "cos4"
## [10] "sin5" "cos5" "sin6" "cos6" "sin7" "cos7" "sin8" "cos8" "sin9"
## [19] "cos9" "sin10" "cos10" "sin11" "cos11"
##
## attr(,"class")
## [1] "basisfd"
##
## $fdnames
## $fdnames[[1]]
## [1] "const" "sin1" "cos1" "sin2" "cos2" "sin3" "cos3" "sin4" "cos4"
## [10] "sin5" "cos5" "sin6" "cos6" "sin7" "cos7" "sin8" "cos8" "sin9"
## [19] "cos9" "sin10" "cos10" "sin11" "cos11"
##
## $fdnames[[2]]
## [1] "reps 1"
##
## $fdnames[[3]]
## [1] "values"
##
##
## attr(,"class")
## [1] "fd"
##
## $Lfd
## $call
## Lfd(nderiv = m, bwtlist = bwtlist)
##
## $nderiv
## [1] 3
##
## $bwtlist
## $bwtlist[[1]]
## $coefs
## reps 1
## time 0
##
## $basis
## $call
## basisfd(type = type, rangeval = rangeval, nbasis = nbasis, params = params,
## dropind = dropind, quadvals = quadvals, values = values,
## basisvalues = basisvalues)
##
## $type
## [1] "const"
##
## $rangeval
## [1] 0 24
##
## $nbasis
## [1] 1
##
## $params
## [1] 0
##
## $dropind
## numeric(0)
##
## $quadvals
## numeric(0)
##
## $values
## list()
##
## $basisvalues
## list()
##
## $names
## [1] "const"
##
## attr(,"class")
## [1] "basisfd"
##
## $fdnames
## $fdnames$args
## [1] "time"
##
## $fdnames$reps
## [1] "reps 1"
##
## $fdnames$funs
## [1] "values"
##
##
## attr(,"class")
## [1] "fd"
##
## $bwtlist[[2]]
## $coefs
## reps 1
## time 0.06853892
##
## $basis
## $call
## basisfd(type = type, rangeval = rangeval, nbasis = nbasis, params = params,
## dropind = dropind, quadvals = quadvals, values = values,
## basisvalues = basisvalues)
##
## $type
## [1] "const"
##
## $rangeval
## [1] 0 24
##
## $nbasis
## [1] 1
##
## $params
## [1] 0
##
## $dropind
## numeric(0)
##
## $quadvals
## numeric(0)
##
## $values
## list()
##
## $basisvalues
## list()
##
## $names
## [1] "const"
##
## attr(,"class")
## [1] "basisfd"
##
## $fdnames
## $fdnames$args
## [1] "time"
##
## $fdnames$reps
## [1] "reps 1"
##
## $fdnames$funs
## [1] "values"
##
##
## attr(,"class")
## [1] "fd"
##
## $bwtlist[[3]]
## $coefs
## reps 1
## time 0
##
## $basis
## $call
## basisfd(type = type, rangeval = rangeval, nbasis = nbasis, params = params,
## dropind = dropind, quadvals = quadvals, values = values,
## basisvalues = basisvalues)
##
## $type
## [1] "const"
##
## $rangeval
## [1] 0 24
##
## $nbasis
## [1] 1
##
## $params
## [1] 0
##
## $dropind
## numeric(0)
##
## $quadvals
## numeric(0)
##
## $values
## list()
##
## $basisvalues
## list()
##
## $names
## [1] "const"
##
## attr(,"class")
## [1] "basisfd"
##
## $fdnames
## $fdnames$args
## [1] "time"
##
## $fdnames$reps
## [1] "reps 1"
##
## $fdnames$funs
## [1] "values"
##
##
## attr(,"class")
## [1] "fd"
##
##
## attr(,"class")
## [1] "Lfd"
##
## $lambda
## [1] 0
##
## $estimate
## [1] TRUE
##
## $penmat
## NULL
##
## attr(,"class")
## [1] "fdPar"
## $fd
## $coefs
## [,1]
## [1,] -135.9295
##
## $basis
## $call
## basisfd(type = type, rangeval = rangeval, nbasis = nbasis, params = params,
## dropind = dropind, quadvals = quadvals, values = values,
## basisvalues = basisvalues)
##
## $type
## [1] "const"
##
## $rangeval
## [1] 0 24
##
## $nbasis
## [1] 1
##
## $params
## [1] 0
##
## $dropind
## numeric(0)
##
## $quadvals
## numeric(0)
##
## $values
## list()
##
## $basisvalues
## list()
##
## $names
## [1] "const"
##
## attr(,"class")
## [1] "basisfd"
##
## $fdnames
## $fdnames[[1]]
## [1] "const"
##
## $fdnames[[2]]
## [1] "reps 1"
##
## $fdnames[[3]]
## [1] "values"
##
##
## attr(,"class")
## [1] "fd"
##
## $Lfd
## $call
## Lfd(nderiv = m, bwtlist = bwtlist)
##
## $nderiv
## [1] 0
##
## $bwtlist
## list()
##
## attr(,"class")
## [1] "Lfd"
##
## $lambda
## [1] 0
##
## $estimate
## [1] TRUE
##
## $penmat
## NULL
##
## attr(,"class")
## [1] "fdPar"# Plot coefficients for model with constant
par(mfrow = c(2, 3)) # Adjust layout if needed
for (i in seq_along(fRegress_no_const$betaestlist)) {
plot(fRegress_no_const$betaestlist[[i]]$fd,
main = paste("Coefficient Function for Predictor", i))
}
par(mfrow = c(1, 1))
# Create a time grid for evaluation (e.g., 24 hours for daily curves)
time_grid <- seq(0, 24, length.out = 100)
# Evaluate the functional objects for predicted and observed data
predicted_with_const_vals <- eval.fd(time_grid, pred_with_const)
predicted_no_const_vals <- eval.fd(time_grid, pred_no_const)
observed_vals <- eval.fd(time_grid, cnt_fd_reg)
# Combine data into a single data frame for ggplot
all_days_data <- data.frame(
time = rep(time_grid, ncol(observed_vals)),
day = factor(rep(1:ncol(observed_vals), each = length(time_grid))),
observed = as.vector(observed_vals),
predicted_with_const = as.vector(predicted_with_const_vals),
predicted_no_const = as.vector(predicted_no_const_vals)
)
# Plot all days
library(ggplot2)
p_all_days <- ggplot(all_days_data, aes(x = time)) +
geom_line(aes(y = observed, color = "Observed")) +
geom_line(aes(y = predicted_with_const, color = "Predicted (With Const)"), linetype = "dashed") +
geom_line(aes(y = predicted_no_const, color = "Predicted (No Const)"), linetype = "dotted") +
labs(
title = "Comparison of Observed and Predicted Curves (All Days)",
x = "Time (Hour of Day)",
y = "Count",
color = "Legend"
) +
theme_minimal()
print(p_all_days)
# Plot for Observed vs Predicted (With Constant)
p_all_days_with_const <- ggplot(all_days_data, aes(x = time)) +
geom_line(aes(y = observed, color = "Observed")) +
geom_line(aes(y = predicted_with_const, color = "Predicted (With Const)"), linetype = "dashed") +
labs(
title = "Comparison of Observed and Predicted Curves (With Constant, All Days)",
x = "Time (Hour of Day)",
y = "Count",
color = "Legend"
) +
theme_minimal()
# Print the plot
print(p_all_days_with_const)
# Plot for Observed vs Predicted (Without Constant)
p_all_days_no_const <- ggplot(all_days_data, aes(x = time)) +
geom_line(aes(y = observed, color = "Observed")) +
geom_line(aes(y = predicted_no_const, color = "Predicted (No Const)"), linetype = "dotted") +
labs(
title = "Comparison of Observed and Predicted Curves (Without Constant, All Days)",
x = "Time (Hour of Day)",
y = "Count",
color = "Legend"
) +
theme_minimal()
# Print the plot
print(p_all_days_no_const)
# Sample 3 random days (or specify the days explicitly)
sampled_days <- sample(1:ncol(observed_vals), 3) # Random selection
# sampled_days <- c(1, 5, 10) # Specify days manually if needed
# Filter data for sampled days
sampled_days_data <- all_days_data %>% filter(day %in% sampled_days)
# Plot sampled days
p_sampled_days <- ggplot(sampled_days_data, aes(x = time)) +
geom_line(aes(y = observed, color = "Observed")) +
geom_line(aes(y = predicted_with_const, color = "Predicted (With Const)"), linetype = "dashed") +
geom_line(aes(y = predicted_no_const, color = "Predicted (No Const)"), linetype = "dotted") +
facet_wrap(~ day, scales = "free_y", ncol = 3) + # Separate panels for each day
labs(
title = "Closer Comparison of Observed and Predicted Curves (Sampled Days)",
x = "Time (Hour of Day)",
y = "Count",
color = "Legend"
) +
theme_minimal()
print(p_sampled_days)
Depth Analysis and Functional Boxplots
# Include t1 in the list of variables
depth_vars <- c("cnt", "t1", "t2", "hum", "wind_speed")
basis_depth <- create.fourier.basis(rangeval = c(0, 24), nbasis = nbasis)
smoothed_fd_list <- lapply(depth_vars, function(var) {
variable_data <- hourly_data %>%
select(matches(paste0("^", var))) %>%
mutate(across(everything(), as.numeric)) %>%
as.matrix()
variable_data[is.na(variable_data)] <- 0
smooth.basis(seq(0, 23, length.out = ncol(variable_data)), t(variable_data), basis_depth)$fd
})
names(smoothed_fd_list) <- depth_vars
# 5 variables, arrange them in a 2x3 grid
# This results in only one blank plot instead of four.
par(mfrow = c(2, 3))
# Plot boxplots for each variable (5 plots)
boxplot(smoothed_fd_list$cnt, method = "MBD", main = "Functional Boxplot: cnt")## $depth
## [1] 0.05972722 0.17190444 0.28877557 0.28417012 0.30944814 0.29283033
## [7] 0.38988070 0.30623366 0.25326248 0.24568375 0.36392420 0.37636331
## [13] 0.36045760 0.41049164 0.32145312 0.25129783 0.34644541 0.35142264
## [19] 0.36221771 0.38386860 0.38585080 0.32051099 0.27494309 0.35637007
## [25] 0.41190082 0.33781975 0.34559842 0.38179704 0.20143582 0.22507467
## [31] 0.31752200 0.29177078 0.34730157 0.37323732 0.38630319 0.29503188
## [37] 0.27104010 0.36687285 0.35316887 0.38719660 0.40617345 0.31424354
## [43] 0.33011871 0.30454884 0.24754119 0.42227798 0.43229664 0.28043887
## [49] 0.36826469 0.32577343 0.20875005 0.28883364 0.36109879 0.39111729
## [55] 0.33901359 0.45252503 0.27954234 0.24573150 0.31208708 0.36386676
## [61] 0.40470968 0.43494827 0.46733579 0.29055553 0.26884134 0.36280313
## [67] 0.41479002 0.45466457 0.46439642 0.47285736 0.34622172 0.20174382
## [73] 0.32544376 0.44018990 0.44119256 0.40050634 0.43606893 0.34717596
## [79] 0.30714128 0.42143366 0.38417704 0.41643331 0.36373640 0.46635565
## [85] 0.34529914 0.12746548 0.34641995 0.34981742 0.37250455 0.42378909
## [91] 0.25412528 0.29278669 0.33421589 0.25830676 0.44170164 0.46803130
## [97] 0.43113304 0.40555593 0.31656339 0.27539277 0.41725136 0.39700089
## [103] 0.37206493 0.46238613 0.46444436 0.31897991 0.32280631 0.39874348
## [109] 0.37352662 0.39257386 0.43051341 0.43759990 0.32488046 0.33447104
## [115] 0.39635098 0.39574552 0.34060905 0.40452897 0.42518907 0.33819729
## [121] 0.35337779 0.23565559 0.42831476 0.38040360 0.44919543 0.44963208
## [127] 0.31204660 0.29675039 0.42546196 0.38672919 0.35057639 0.21538649
## [133] 0.42005977 0.26423533 0.31658302 0.36152639 0.40416238 0.38430507
## [139] 0.38051774 0.37801091 0.34139823 0.24979013 0.34739100 0.40654883
## [145] 0.37763311 0.38007775 0.38661170 0.30753741 0.31763797 0.40739627
## [151] 0.42995538 0.36405901 0.29365446 0.32921294 0.27281884 0.26947496
## [157] 0.43518554 0.39022837 0.37047968 0.29919273 0.34119009 0.26196376
## [163] 0.35064129 0.40539392 0.32729174 0.30417972 0.28560647 0.32364317
## [169] 0.26842491 0.29282859 0.39326537 0.32009979 0.29003324 0.28073507
## [175] 0.30597037 0.19946757 0.32818148 0.30968715 0.21332673 0.21592222
## [181] 0.27593892 0.26021126 0.18167533 0.27913942 0.32569153 0.31737238
## [187] 0.25306721 0.09278057 0.20371148 0.15961339 0.28347761 0.38967364
## [193] 0.39201880 0.39291358 0.27553396 0.26213854 0.19800904 0.20017778
## [199] 0.33940274 0.26366365 0.31341039 0.27250442 0.26492142 0.25888528
## [205] 0.21806651 0.44557524 0.36594392 0.35331686 0.35030951 0.30419204
## [211] 0.16725818 0.18368751 0.40827350 0.33249107 0.27756140 0.11497513
## [217] 0.25555416 0.17017398 0.20625264 0.39949258 0.40240319 0.34234860
## [223] 0.36003431 0.40696578 0.23718177 0.30094806 0.40511639 0.40993955
## [229] 0.39539849 0.42212509 0.33016977 0.13968862 0.31136860 0.29436304
## [235] 0.35335081 0.29763486 0.42886088 0.38585511 0.26283506 0.29619915
## [241] 0.20966765 0.41037516 0.43726244 0.42646013 0.44200396 0.34682810
## [247] 0.31426811 0.43668619 0.41541365 0.39745435 0.36159467 0.34470233
## [253] 0.28638280 0.30744583 0.39417704 0.42905187 0.32926790 0.42068107
## [259] 0.40395836 0.23870684 0.28679467 0.37745518 0.36578417 0.40313652
## [265] 0.45451955 0.36019175 0.25602250 0.29348859 0.43405771 0.40282103
## [271] 0.38819307 0.39971482 0.35595142 0.26346467 0.27889818 0.26987383
## [277] 0.38857639 0.40235351 0.40302078 0.29354626 0.26761072 0.28275193
## [283] 0.42587527 0.25675312 0.41659707 0.43449251 0.44320934 0.30144241
## [289] 0.30076018 0.45386021 0.40980748 0.34093965 0.45004862 0.41794903
## [295] 0.28017767 0.24632972 0.46819822 0.44113144 0.46412835 0.44089413
## [301] 0.44377654 0.26876897 0.23072448 0.44978475 0.41450651 0.42577979
## [307] 0.39659114 0.40571223 0.27264391 0.27808600 0.45056773 0.40482910
## [313] 0.40537811 0.40937339 0.42065435 0.20325545 0.30162283 0.39640508
## [319] 0.28266517 0.34736973 0.36152579 0.45809595 0.28311076 0.27345650
## [325] 0.33675872 0.33079103 0.40980614 0.45600406 0.42849911 0.34405457
## [331] 0.21014267 0.29951342 0.39328626 0.43382623 0.41716783 0.46271246
## [337] 0.32597827 0.26983442 0.40642177 0.40697776 0.46566956 0.37092969
## [343] 0.40931067 0.32134425 0.25682133 0.37539980 0.35269652 0.43382315
## [349] 0.43092353 0.42079644 0.30740913 0.23774878 0.34456433 0.33195289
## [355] 0.39788373 0.22756877 0.34523732 0.23429977 0.21051736 0.30278925
## [361] 0.35615704 0.27262773 0.28311121 0.11561773 0.13959485 0.10260541
## [367] 0.28994963 0.32796941 0.32765956 0.28922521 0.35013250 0.27206610
## [373] 0.25341236 0.28223946 0.32291566 0.34113012 0.31690571 0.39089060
## [379] 0.33096570 0.24558630 0.30851762 0.33900108 0.36822273 0.39647295
## [385] 0.31907984 0.29743990 0.30401931 0.42201885 0.29145203 0.31909695
## [391] 0.42249577 0.40945016 0.33036830 0.22757453 0.35136739 0.38875269
## [397] 0.38551705 0.44147191 0.41693627 0.28798198 0.24423816 0.23998191
## [403] 0.33745386 0.38876627 0.42450680 0.42124530 0.19797126 0.26521706
## [409] 0.31657939 0.39362494 0.26676469 0.39871721 0.41927649 0.24005044
## [415] 0.28715884 0.29156958 0.39772706 0.42010356 0.40916254 0.44420681
## [421] 0.32772902 0.26399339 0.32801275 0.28504279 0.27987705 0.39379156
## [427] 0.42299220 0.27545288 0.25386541 0.30864271 0.32620548 0.29908516
## [433] 0.40647677 0.37929173 0.32437427 0.24878476 0.39399829 0.38609653
## [439] 0.40372469 0.46207613 0.42812721 0.36374504 0.31172721 0.43266690
## [445] 0.44104594 0.45507376 0.32681490 0.31238513 0.23385637 0.11673880
## [451] 0.17966515 0.32712884 0.43704220 0.47976358 0.46803056 0.22900493
## [457] 0.25712038 0.42806187 0.42870064 0.40024098 0.27812084 0.46116765
## [463] 0.33340359 0.25774183 0.37836314 0.40921304 0.43876283 0.43855907
## [469] 0.29612972 0.38088459 0.30054443 0.40192057 0.40182197 0.39889737
## [475] 0.44971680 0.35561663 0.35679604 0.27939765 0.35506636 0.37803069
## [481] 0.36918950 0.43337827 0.45201973 0.35345427 0.32765614 0.33255255
## [487] 0.43166873 0.32743238 0.28658044 0.26606038 0.17771152 0.22209081
## [493] 0.37614489 0.31628170 0.33964298 0.30221593 0.38406668 0.31746430
## [499] 0.30599260 0.35720757 0.33001841 0.35725614 0.38502782 0.35820959
## [505] 0.28834563 0.29893988 0.38603761 0.34631661 0.42059082 0.33926221
## [511] 0.35310842 0.27049843 0.27777009 0.37027855 0.31756197 0.45029272
## [517] 0.47487856 0.45702353 0.27971511 0.22641746 0.33174531 0.30550990
## [523] 0.33177206 0.23866268 0.25299623 0.27693527 0.31940427 0.38794760
## [529] 0.32680793 0.32131887 0.42177038 0.41221731 0.28601499 0.23571696
## [535] 0.33081490 0.30478284 0.33111054 0.36098610 0.15939835 0.29206108
## [541] 0.29491395 0.40566199 0.35490557 0.32419767 0.36357238 0.39814757
## [547] 0.26621108 0.28362846 0.34715313 0.26146611 0.22523965 0.24985269
## [553] 0.30527092 0.23641912 0.30688587 0.41502766 0.29437217 0.33979022
## [559] 0.25869929 0.30502508 0.18109444 0.20353050 0.24054009 0.17844507
## [565] 0.20021063 0.22881215 0.31127052 0.15089378 0.22979748 0.29441859
## [571] 0.23349504 0.33300179 0.24482536 0.33303852 0.16052863 0.16841747
## [577] 0.35961577 0.36615087 0.28695995 0.26936055 0.18447636 0.15294982
## [583] 0.23573563 0.35773713 0.30086861 0.30681113 0.33352174 0.26014361
## [589] 0.14181588 0.16954014 0.25634085 0.25735349 0.25091127 0.29650403
## [595] 0.38001664 0.30511199 0.32030904 0.39152210 0.23071903 0.27349825
## [601] 0.32959267 0.28463876 0.24975287 0.28823379 0.27643053 0.31460572
## [607] 0.31546644 0.23182911 0.15711016 0.32902955 0.43652945 0.34796673
## [613] 0.29289241 0.26400456 0.31046772 0.27233881 0.26753398 0.32230008
## [619] 0.22710708 0.21268714 0.21936826 0.35163802 0.31927547 0.31593322
## [625] 0.41291828 0.33382276 0.27814659 0.30147936 0.28592070 0.23404198
## [631] 0.30031614 0.38116120 0.34380194 0.31684690 0.36659008 0.34428190
## [637] 0.28437410 0.31977101 0.38632137 0.35462748 0.34912057 0.38356435
## [643] 0.38569450 0.30585011 0.26434477 0.40383434 0.35978216 0.41981029
## [649] 0.44082352 0.40635168 0.25961074 0.28265945 0.42522176 0.37494905
## [655] 0.38185378 0.41691430 0.43262059 0.31878954 0.28592942 0.44475445
## [661] 0.41313162 0.40218407 0.41111383 0.36033722 0.27599781 0.27299198
## [667] 0.39988612 0.38951968 0.40018194 0.43453808 0.36894581 0.37301260
## [673] 0.25103547 0.22526525 0.35533105 0.32107866 0.40325245 0.47360987
## [679] 0.19261553 0.28820359 0.37220832 0.41682851 0.40750926 0.43493217
## [685] 0.46961749 0.28542720 0.18418974 0.26252970 0.39739557 0.43856787
## [691] 0.43359683 0.47001532 0.34078862 0.29887409 0.41169049 0.40222541
## [697] 0.42199588 0.44448753 0.45762921 0.34842557 0.28285483 0.41937795
## [703] 0.42717743 0.45580642 0.43705931 0.43575489 0.27560866 0.30010218
## [709] 0.41153337 0.41424597 0.46970132 0.44725991 0.43500493 0.33818823
## [715] 0.30266679 0.43360429 0.46414316 0.45660168 0.44242740 0.33495322
## [721] 0.16320030 0.25003422 0.24333555 0.21351469 0.26631684 0.31129482
## [727] 0.25732763 0.20387638 0.07427358 0.19690292 0.26524074
##
## $outpoint
## [1] 188
##
## $medcurve
## [1] 454## $depth
## [1] 0.018804589 0.064633132 0.373194836 0.328668444 0.318960981 0.408865489
## [7] 0.479287109 0.350525849 0.285265999 0.485407826 0.336902057 0.268359738
## [13] 0.371261654 0.205032634 0.113976690 0.097089818 0.048168244 0.050250637
## [19] 0.079561700 0.057731587 0.106644145 0.195457490 0.225226401 0.368450077
## [25] 0.298797355 0.268117942 0.104076080 0.111068354 0.097392434 0.101739088
## [31] 0.021225409 0.044578548 0.067042856 0.068061250 0.059488172 0.110596414
## [37] 0.190693027 0.216265210 0.191642372 0.172091752 0.115791952 0.311872124
## [43] 0.313423984 0.265126512 0.247633071 0.216712910 0.267272880 0.318517985
## [49] 0.256840511 0.151685625 0.200889369 0.280737998 0.311386268 0.422453537
## [55] 0.390980837 0.235910436 0.403333112 0.396979060 0.236243064 0.269713960
## [61] 0.284684330 0.329671475 0.418515221 0.469804235 0.446862698 0.377958177
## [67] 0.393335938 0.389112418 0.447175556 0.367207662 0.247508007 0.248960594
## [73] 0.297883059 0.379174691 0.294919460 0.308322598 0.330691687 0.324944677
## [79] 0.209381960 0.364653189 0.319403087 0.251210911 0.381156818 0.397269171
## [85] 0.490576563 0.459775314 0.432837072 0.371603020 0.393377462 0.405247601
## [91] 0.434373513 0.368577479 0.397009637 0.394187280 0.429012259 0.474717584
## [97] 0.443721212 0.431356858 0.478172791 0.450326774 0.476865242 0.368118072
## [103] 0.314642854 0.488930653 0.463301268 0.425637698 0.380713590 0.427505857
## [109] 0.449955480 0.456485007 0.453661351 0.454830757 0.474454544 0.414347233
## [115] 0.409145423 0.433660732 0.445500433 0.458451289 0.415712661 0.452773198
## [121] 0.447300712 0.458643110 0.437936166 0.498547089 0.493238279 0.480786647
## [127] 0.456173141 0.438388216 0.324455305 0.464392263 0.450657972 0.441914708
## [133] 0.486438967 0.436710440 0.478867232 0.495286718 0.476686576 0.481917951
## [139] 0.406207921 0.371160043 0.440095976 0.401700848 0.478876462 0.449764030
## [145] 0.428541024 0.445258071 0.492837168 0.479909115 0.474510809 0.490538119
## [151] 0.405424449 0.393950141 0.338866595 0.200276529 0.425292119 0.431780085
## [157] 0.476809506 0.496498093 0.471102442 0.329719279 0.233352735 0.307596523
## [163] 0.447247379 0.468601887 0.300101027 0.182914776 0.306683266 0.372875373
## [169] 0.344200918 0.327766671 0.452426775 0.397874774 0.245331035 0.181944650
## [175] 0.137848211 0.185431779 0.209063591 0.131322227 0.072231418 0.006460998
## [181] 0.076014289 0.161697067 0.061550602 0.242073196 0.221485249 0.195665842
## [187] 0.294822849 0.291507206 0.117647346 0.092320800 0.242971731 0.219199234
## [193] 0.190070910 0.171095929 0.152361316 0.138212909 0.243395255 0.163737603
## [199] 0.167496539 0.140656883 0.197694039 0.263015721 0.376112143 0.350984273
## [205] 0.393885796 0.290978972 0.343297197 0.399906191 0.404249115 0.338831649
## [211] 0.301878643 0.178895443 0.172478019 0.269877542 0.249379776 0.194581145
## [217] 0.197009961 0.178705459 0.108320576 0.131782748 0.264082736 0.246507732
## [223] 0.200980247 0.219754278 0.313125375 0.350563440 0.281344455 0.291704138
## [229] 0.287755904 0.188933585 0.059410978 0.048933709 0.217104919 0.395239647
## [235] 0.365315744 0.261867041 0.352425059 0.341954701 0.253976909 0.289411657
## [241] 0.403322545 0.418185571 0.453276875 0.482524037 0.463067079 0.491292823
## [247] 0.432390402 0.466433087 0.459966523 0.412557186 0.387635173 0.311816509
## [253] 0.278136370 0.476245908 0.451901594 0.468925979 0.466903868 0.447257435
## [259] 0.441732925 0.408685552 0.417553305 0.439860043 0.499194848 0.450561677
## [265] 0.412995692 0.485080366 0.479911137 0.472326396 0.456273110 0.453690073
## [271] 0.467090985 0.467274002 0.476001134 0.494554038 0.475722007 0.383370565
## [277] 0.354202972 0.433777835 0.489151075 0.462932284 0.490005071 0.493297253
## [283] 0.490278001 0.358635809 0.469903499 0.486377293 0.494064249 0.493276583
## [289] 0.494675836 0.493920668 0.489735137 0.470989967 0.466137753 0.496443377
## [295] 0.475649220 0.358846035 0.478247601 0.462412984 0.474784954 0.476400621
## [301] 0.442734110 0.492284805 0.421505565 0.472031322 0.496532149 0.485247686
## [307] 0.414696754 0.380841734 0.389444127 0.455712370 0.468106856 0.385057229
## [313] 0.432444941 0.462580722 0.435353899 0.442728646 0.402114659 0.460076938
## [319] 0.380386900 0.388366981 0.463854104 0.388087919 0.130257495 0.096508891
## [325] 0.189703289 0.408760362 0.407649708 0.435875972 0.467486965 0.360544546
## [331] 0.498057661 0.465346052 0.484563617 0.479562795 0.479045796 0.477140649
## [337] 0.489778109 0.475316917 0.461745345 0.454774278 0.430899001 0.480810025
## [343] 0.431110628 0.394894662 0.373247808 0.423692056 0.481950764 0.447838677
## [349] 0.454508343 0.481750900 0.410495832 0.434991640 0.480028659 0.452849121
## [355] 0.480119369 0.429381243 0.445768279 0.431208760 0.451194583 0.487972281
## [361] 0.468528728 0.466821747 0.423647638 0.269852902 0.469102243 0.390043672
## [367] 0.411086843 0.401847119 0.301190585 0.337425921 0.329988944 0.435056394
## [373] 0.311124870 0.293032867 0.244110086 0.219873877 0.137043782 0.111389701
## [379] 0.074125859 0.076055033 0.080213735 0.022771220 0.030408877 0.129681540
## [385] 0.364538321 0.402234500 0.498083377 0.481755408 0.458423208 0.473964801
## [391] 0.353504431 0.486817998 0.324552166 0.425704056 0.488618704 0.383853536
## [397] 0.247447456 0.460310821 0.471008623 0.474026400 0.338943938 0.306005923
## [403] 0.245607165 0.264105483 0.189724886 0.149729436 0.131782534 0.130431012
## [409] 0.097882419 0.109221337 0.160650861 0.178321019 0.269536073 0.451198043
## [415] 0.491772027 0.379549436 0.247307689 0.110113267 0.076286105 0.123802397
## [421] 0.156584029 0.169939039 0.211747850 0.424869227 0.245737442 0.301122909
## [427] 0.221527134 0.106735644 0.128177587 0.156724038 0.222853228 0.353544860
## [433] 0.277652889 0.192985201 0.291452545 0.196982789 0.272274259 0.267192466
## [439] 0.304966080 0.263535095 0.216990858 0.308979903 0.355970151 0.403882098
## [445] 0.384097160 0.409050835 0.416588486 0.469122967 0.377949595 0.264767278
## [451] 0.388164139 0.391336332 0.423513789 0.386743253 0.387617037 0.463571015
## [457] 0.484102086 0.480770801 0.488754854 0.424773673 0.334391891 0.432705477
## [463] 0.394281154 0.386842406 0.496763258 0.485705424 0.464584984 0.487793059
## [469] 0.467798692 0.269507982 0.340083341 0.432537006 0.472762230 0.411271716
## [475] 0.468665657 0.404454963 0.334483325 0.319646071 0.339940800 0.253864741
## [481] 0.270568865 0.359638961 0.385999459 0.406106059 0.450350866 0.494765201
## [487] 0.487067800 0.466037116 0.436255773 0.405453022 0.311989562 0.223618202
## [493] 0.298432945 0.394142407 0.387217819 0.312921569 0.469711297 0.458076340
## [499] 0.474814158 0.462405544 0.459888559 0.481213083 0.464817058 0.395886821
## [505] 0.389545487 0.431563550 0.442059838 0.489153867 0.482172345 0.456882434
## [511] 0.430795321 0.459297650 0.468743074 0.456094110 0.484229784 0.487522559
## [517] 0.483573611 0.488958252 0.424319229 0.378949148 0.395997032 0.302729169
## [523] 0.264899457 0.407788019 0.299676835 0.263515752 0.320351013 0.353517669
## [529] 0.383430949 0.374209723 0.390278131 0.419994284 0.434718794 0.338437895
## [535] 0.263145395 0.237366297 0.193152596 0.163973332 0.150761857 0.381847511
## [541] 0.362408392 0.268290299 0.358823993 0.410978488 0.320014081 0.355162308
## [547] 0.415850156 0.350415787 0.306103183 0.280589519 0.290168940 0.189307569
## [553] 0.176597646 0.138338659 0.199090435 0.275554538 0.403489810 0.409831542
## [559] 0.325111552 0.305375420 0.072716086 0.033464170 0.032562991 0.032146277
## [565] 0.022592136 0.114835668 0.097760880 0.063189962 0.107791675 0.129542125
## [571] 0.167937884 0.124167838 0.185350921 0.156793477 0.188549154 0.226415158
## [577] 0.357342073 0.196094858 0.140775592 0.197102639 0.183540084 0.137294549
## [583] 0.092257336 0.255448597 0.368102050 0.367926074 0.210766239 0.129916397
## [589] 0.136855635 0.213278084 0.265716613 0.275521299 0.210871014 0.191559166
## [595] 0.279701876 0.275469579 0.224540625 0.122779699 0.064845724 0.010572201
## [601] 0.021263992 0.038570649 0.087894141 0.167036418 0.195269194 0.133023039
## [607] 0.167844576 0.190533286 0.170332319 0.230222214 0.204577911 0.063905248
## [613] 0.063083341 0.156556310 0.161837020 0.222469531 0.351875015 0.122641026
## [619] 0.025148400 0.031768999 0.058191393 0.320076673 0.412458849 0.314185469
## [625] 0.276986474 0.340240958 0.280286764 0.317831375 0.415075720 0.200430825
## [631] 0.319957621 0.429535529 0.316764202 0.280019511 0.306780554 0.432390773
## [637] 0.482344295 0.491890597 0.468080527 0.438740483 0.464968403 0.482150934
## [643] 0.466064344 0.477279247 0.494954601 0.462781207 0.466570721 0.498245280
## [649] 0.488410602 0.496812908 0.490672905 0.460208180 0.469318091 0.492930291
## [655] 0.498851320 0.494945463 0.482828648 0.466888189 0.466167393 0.497119171
## [661] 0.495992922 0.487674007 0.480610597 0.451924387 0.463215789 0.432616919
## [667] 0.490953702 0.463834937 0.355527201 0.355141871 0.422340692 0.298342568
## [673] 0.214051491 0.271121367 0.156886952 0.340209944 0.366397454 0.319118320
## [679] 0.408164138 0.388717476 0.468060582 0.464920822 0.459416136 0.453564907
## [685] 0.235484742 0.223931431 0.319163072 0.450573923 0.456400484 0.416529076
## [691] 0.440098852 0.367025656 0.254424173 0.341281078 0.263567230 0.075668367
## [697] 0.070966404 0.149311450 0.289170501 0.264929458 0.144077348 0.104329640
## [703] 0.313939944 0.496481394 0.494621352 0.482805929 0.485731501 0.356539175
## [709] 0.365684292 0.489587549 0.442967547 0.419801785 0.434477054 0.318723425
## [715] 0.363521792 0.305049611 0.239293774 0.467888809 0.369599250 0.452666419
## [721] 0.413124429 0.496448396 0.329725922 0.116121212 0.047246507 0.109152158
## [727] 0.096735204 0.329591052 0.287822699 0.097364715 0.092627981
##
## $outpoint
## integer(0)
##
## $medcurve
## [1] 263## $depth
## [1] 0.080880530 0.122025997 0.407726420 0.311927035 0.304977844 0.391201361
## [7] 0.476541094 0.308699820 0.203915236 0.471070557 0.302749421 0.208510996
## [13] 0.309150387 0.179290653 0.106577944 0.155008616 0.102847172 0.060519368
## [19] 0.072513476 0.115441392 0.096984088 0.169970961 0.199392578 0.346647276
## [25] 0.289327218 0.201685560 0.031096415 0.077533048 0.038327498 0.033989435
## [31] 0.012016261 0.035149804 0.047870303 0.037441562 0.017755544 0.124119746
## [37] 0.217801484 0.273837975 0.280002470 0.225047614 0.172818551 0.293068927
## [43] 0.334279555 0.323840210 0.281896215 0.247456288 0.264466257 0.276445002
## [49] 0.300887458 0.137669971 0.150717809 0.199390166 0.232646345 0.423380553
## [55] 0.383555150 0.230229821 0.368195896 0.345377178 0.141722961 0.204175734
## [61] 0.245228616 0.314556216 0.410047669 0.463740145 0.458189278 0.383594317
## [67] 0.399158500 0.403166877 0.452530149 0.358039343 0.222184049 0.234040755
## [73] 0.305708873 0.404547325 0.309499870 0.300145148 0.351341772 0.313056855
## [79] 0.235038146 0.372067604 0.362294637 0.298787512 0.354471847 0.389144071
## [85] 0.491297211 0.455045946 0.402735709 0.359220530 0.320951162 0.405984484
## [91] 0.453381074 0.388585669 0.438976528 0.402573445 0.463886443 0.481145296
## [97] 0.455004218 0.438203149 0.474771196 0.456119242 0.485397315 0.375372922
## [103] 0.327029668 0.483257256 0.444686950 0.418430392 0.388911626 0.451037495
## [109] 0.467007232 0.460489904 0.458927840 0.461635254 0.474547472 0.428214677
## [115] 0.416885842 0.410001098 0.429112089 0.446983531 0.407846158 0.444488310
## [121] 0.447287205 0.458636987 0.438005335 0.497232758 0.493087341 0.481212350
## [127] 0.456236410 0.438059810 0.322344004 0.464490823 0.449795061 0.423452292
## [133] 0.485114264 0.436776780 0.478359324 0.491050368 0.457990453 0.472915718
## [139] 0.404763118 0.370414031 0.440103305 0.401715617 0.478978250 0.449855232
## [145] 0.428000730 0.445344347 0.492143637 0.473315085 0.473408375 0.487519053
## [151] 0.405620722 0.394046724 0.346171885 0.211339893 0.424766474 0.431977472
## [157] 0.476055293 0.493318970 0.470665114 0.334896143 0.231173480 0.307504124
## [163] 0.447328265 0.468672216 0.305005452 0.180194680 0.306715336 0.372964219
## [169] 0.344240465 0.327674792 0.452421088 0.397898217 0.251672767 0.182689437
## [175] 0.135538456 0.197142938 0.206092236 0.135425387 0.075343709 0.006407414
## [181] 0.073723942 0.160462853 0.060243581 0.242467767 0.231195355 0.197478878
## [187] 0.294660306 0.301332143 0.132615035 0.093448366 0.240684427 0.224698687
## [193] 0.186656921 0.167446276 0.149389999 0.135256286 0.253424184 0.173058631
## [199] 0.171010776 0.138864230 0.200868774 0.263155526 0.376176961 0.351038590
## [205] 0.393942655 0.290893132 0.343469824 0.399900161 0.404312079 0.340413566
## [211] 0.304606410 0.184771645 0.170560237 0.269661508 0.248189516 0.191574464
## [217] 0.194117860 0.176020319 0.112126825 0.127529309 0.262911280 0.244128362
## [223] 0.199679460 0.218237151 0.313059694 0.350565100 0.280344504 0.291707997
## [229] 0.286665075 0.185456521 0.058534977 0.048969741 0.216809326 0.395341759
## [235] 0.365243077 0.261827521 0.352515677 0.341024810 0.254018479 0.289475186
## [241] 0.403447414 0.418236112 0.453338465 0.482582779 0.463092479 0.491771285
## [247] 0.434376992 0.466567010 0.460041778 0.412742800 0.387739707 0.324838436
## [253] 0.291994824 0.476503178 0.452089778 0.469052128 0.467027373 0.447403901
## [259] 0.441856717 0.415677335 0.417704150 0.439970913 0.499186619 0.450143246
## [265] 0.413096811 0.485205346 0.482624637 0.473384013 0.456941946 0.453777156
## [271] 0.467099808 0.467323031 0.476877849 0.494884291 0.477215542 0.383443501
## [277] 0.354129433 0.433795907 0.488832010 0.477554200 0.491362855 0.493838288
## [283] 0.492338538 0.370099012 0.477416343 0.484856289 0.494123473 0.493395069
## [289] 0.494771204 0.494901880 0.495612472 0.470302614 0.466229967 0.496481710
## [295] 0.475658645 0.386826959 0.479098415 0.462466420 0.474871416 0.476636220
## [301] 0.442852206 0.493745824 0.452306155 0.478816987 0.496412624 0.485343759
## [307] 0.414773864 0.380880289 0.389477719 0.455420339 0.468270242 0.385158237
## [313] 0.432491326 0.462795902 0.404720518 0.434806740 0.402210398 0.460191825
## [319] 0.392406529 0.409672404 0.462413142 0.340890516 0.056465794 0.103561354
## [325] 0.215119350 0.381927664 0.409007725 0.451621502 0.448656697 0.316597531
## [331] 0.498034506 0.471126451 0.485611336 0.479572554 0.479180721 0.464343132
## [337] 0.489763729 0.475353496 0.461884973 0.447035631 0.429310293 0.480251251
## [343] 0.423357898 0.420955336 0.401179248 0.436501169 0.482061494 0.447905212
## [349] 0.454602412 0.481781848 0.410578509 0.432976830 0.469962037 0.452908763
## [355] 0.478697722 0.421307232 0.448941873 0.431303312 0.451310221 0.488098718
## [361] 0.468233320 0.465020049 0.399585392 0.267531848 0.462785326 0.386700597
## [367] 0.415278247 0.434020643 0.372214988 0.278146630 0.344108936 0.422760347
## [373] 0.286900387 0.278768042 0.149218411 0.194224574 0.065338445 0.053117700
## [379] 0.095826362 0.101875052 0.083174180 0.110013363 0.086260412 0.141860845
## [385] 0.350947722 0.405150266 0.496750233 0.477851509 0.450611355 0.462920149
## [391] 0.316747337 0.486533305 0.271342605 0.421942485 0.482485751 0.318069052
## [397] 0.197661820 0.441851893 0.467675836 0.456974954 0.286306175 0.241807067
## [403] 0.190776334 0.239329658 0.238589082 0.158686768 0.102549566 0.113814240
## [409] 0.078299441 0.147328432 0.163762039 0.188822900 0.230177694 0.453739881
## [415] 0.491837653 0.379284827 0.299215405 0.233956149 0.117637327 0.138560463
## [421] 0.134450818 0.123508205 0.265707253 0.391493494 0.158727883 0.259268270
## [427] 0.228317819 0.089053063 0.121310516 0.149936833 0.212783397 0.336092479
## [433] 0.344620165 0.261749056 0.365386890 0.201089011 0.237693229 0.290637095
## [439] 0.251823008 0.293292762 0.231786208 0.292283135 0.367641947 0.435360309
## [445] 0.429864010 0.445346917 0.400105016 0.462706332 0.380284036 0.281010436
## [451] 0.329520677 0.351016158 0.414600097 0.423005446 0.422868137 0.461925384
## [457] 0.488683792 0.482176604 0.486956893 0.388867541 0.307415909 0.437224220
## [463] 0.391065564 0.384039410 0.493794658 0.488421957 0.476600773 0.489757504
## [469] 0.466915835 0.278198062 0.365420092 0.425445728 0.466940354 0.412284738
## [475] 0.466948926 0.391146905 0.346772702 0.325327353 0.325390725 0.191534424
## [481] 0.280119508 0.351296844 0.351985263 0.424327708 0.466844430 0.492042267
## [487] 0.485699625 0.470123605 0.439858401 0.406297091 0.312127975 0.223347174
## [493] 0.299377039 0.394200110 0.387273740 0.314685009 0.464274482 0.454832250
## [499] 0.485345253 0.462407909 0.460043300 0.481298441 0.464792122 0.395933354
## [505] 0.389613432 0.431749295 0.442169418 0.490731517 0.480627713 0.457075609
## [511] 0.430850194 0.459129485 0.468873462 0.456225621 0.484431503 0.487566839
## [517] 0.482353795 0.487308733 0.424460574 0.376811008 0.398638802 0.300781700
## [523] 0.261541213 0.407987800 0.300926032 0.262164054 0.320418178 0.353606135
## [529] 0.383545019 0.374336113 0.390393834 0.420092397 0.434812492 0.338436633
## [535] 0.259792458 0.235434822 0.188753100 0.166227944 0.192340041 0.381863764
## [541] 0.362442763 0.267424884 0.358292383 0.410922288 0.320036160 0.355186911
## [547] 0.416123568 0.350402688 0.305612976 0.280579658 0.294000460 0.186961819
## [553] 0.172709008 0.134020550 0.195713563 0.275421877 0.403585530 0.409972256
## [559] 0.325243638 0.305356160 0.071377450 0.032451973 0.033633235 0.031411574
## [565] 0.022092995 0.112763145 0.095089990 0.062383863 0.105544669 0.141249388
## [571] 0.164557209 0.121996978 0.186903328 0.153206648 0.199186164 0.224830495
## [577] 0.357379627 0.193461224 0.137312222 0.195888750 0.194930295 0.135392909
## [583] 0.090783747 0.255414754 0.368143760 0.368001812 0.206949924 0.129364081
## [589] 0.134518642 0.220761734 0.271599820 0.277393215 0.209195584 0.194576266
## [595] 0.279628689 0.275525130 0.222731737 0.120404151 0.065661285 0.010464736
## [601] 0.018497743 0.042860305 0.085398326 0.163856404 0.195247718 0.131887606
## [607] 0.164903984 0.201353851 0.206536550 0.230290994 0.200960542 0.061257114
## [613] 0.061940756 0.160086577 0.157391946 0.222063309 0.355446704 0.121110504
## [619] 0.025408898 0.031118087 0.055707862 0.320015667 0.412446075 0.314055294
## [625] 0.276934819 0.340373990 0.281497925 0.317907929 0.415230693 0.206387746
## [631] 0.320108651 0.429698165 0.316914638 0.276717561 0.306913104 0.432384502
## [637] 0.482488887 0.492312506 0.473135547 0.438953855 0.464985807 0.482248027
## [643] 0.466136872 0.477319231 0.495871959 0.478606770 0.484375878 0.497953360
## [649] 0.487824758 0.498048505 0.490686904 0.460066733 0.469422402 0.489793165
## [655] 0.498806725 0.493360049 0.484060729 0.483398387 0.476746431 0.496867745
## [661] 0.497080012 0.487857979 0.480607730 0.452015748 0.463493729 0.454428700
## [667] 0.491820435 0.476218736 0.384239349 0.382273465 0.440195454 0.275919589
## [673] 0.169604704 0.257382168 0.165853329 0.312606410 0.346247465 0.346251927
## [679] 0.417511039 0.426743741 0.474362386 0.465057871 0.459038756 0.427768489
## [685] 0.207678727 0.193758979 0.300236601 0.449467508 0.458495680 0.437972671
## [691] 0.424768896 0.366049890 0.329885375 0.370647998 0.253205404 0.137324643
## [697] 0.075803320 0.155414801 0.366423040 0.276358615 0.158850154 0.168218495
## [703] 0.372153064 0.498235808 0.494700160 0.482975383 0.483730949 0.348323142
## [709] 0.389995024 0.493960244 0.453812261 0.427465650 0.449603055 0.387786249
## [715] 0.414942966 0.378406943 0.263621863 0.463834492 0.373811498 0.450125184
## [721] 0.372084637 0.496586670 0.303258330 0.188190755 0.070244206 0.185760260
## [727] 0.124941615 0.333868175 0.280457350 0.106166749 0.055284792
##
## $outpoint
## integer(0)
##
## $medcurve
## [1] 263## $depth
## [1] 0.02324561 0.06368554 0.38033621 0.38337985 0.38428462 0.38843226
## [7] 0.44756428 0.40720100 0.38288468 0.25870631 0.43785386 0.40383401
## [13] 0.36457619 0.32453866 0.43961398 0.26753398 0.36373959 0.29758692
## [19] 0.30858738 0.46165851 0.38698797 0.46011404 0.44072355 0.44411074
## [25] 0.44887630 0.40254723 0.44511038 0.41030792 0.31188216 0.37613465
## [31] 0.36475657 0.24309144 0.38841059 0.39423922 0.41941846 0.33256907
## [37] 0.45680822 0.33360979 0.32308191 0.44152823 0.41676435 0.30033710
## [43] 0.17886829 0.20622992 0.17095373 0.46114630 0.46213806 0.32450341
## [49] 0.36388950 0.43147536 0.35557707 0.42605959 0.42735885 0.40040028
## [55] 0.23054566 0.41937089 0.26073374 0.35752410 0.34814180 0.37366017
## [61] 0.39891767 0.40259957 0.34097267 0.33553278 0.42452543 0.30742583
## [67] 0.38390935 0.42074674 0.19152793 0.28547076 0.38312695 0.35153389
## [73] 0.32111566 0.22861259 0.36434874 0.42657465 0.39818768 0.45719226
## [79] 0.40817075 0.39037725 0.47178248 0.40827247 0.32253074 0.27794812
## [85] 0.33952256 0.29671447 0.26600992 0.15949413 0.23679748 0.32958776
## [91] 0.15873672 0.36039504 0.46780465 0.45176424 0.42314289 0.43625035
## [97] 0.28728909 0.26954100 0.23260794 0.28285190 0.44234829 0.28336636
## [103] 0.16385094 0.38643581 0.35898749 0.37700056 0.43418603 0.31571080
## [109] 0.31861146 0.41986447 0.45292253 0.40975230 0.36598923 0.29990955
## [115] 0.17913784 0.18605545 0.34610345 0.26189496 0.38657003 0.40189964
## [121] 0.31738526 0.29344658 0.35744105 0.42847410 0.35385196 0.40891343
## [127] 0.25496177 0.40908116 0.35702058 0.17264789 0.23525423 0.23888432
## [133] 0.43484925 0.25333695 0.24831237 0.30386579 0.33769791 0.31295191
## [139] 0.16556393 0.41733737 0.39695279 0.27948675 0.29668619 0.29010370
## [145] 0.19601667 0.09103827 0.24216226 0.27025200 0.28195440 0.33539988
## [151] 0.42245531 0.13314922 0.24280153 0.26534887 0.10961212 0.21776363
## [157] 0.34667678 0.24760100 0.25890643 0.35312353 0.41432049 0.40687842
## [163] 0.39278597 0.36943317 0.25616459 0.41700488 0.19913928 0.26058679
## [169] 0.27599955 0.34147095 0.40476349 0.31133189 0.20677556 0.21549043
## [175] 0.20616992 0.25509372 0.31384273 0.16378668 0.15968898 0.10322193
## [181] 0.39753825 0.32277269 0.24771768 0.28881160 0.25827708 0.23566349
## [187] 0.29317606 0.09774469 0.08800190 0.19768946 0.44315487 0.30245271
## [193] 0.39381308 0.28780612 0.42204735 0.34124456 0.15629719 0.15481165
## [199] 0.30449377 0.25923510 0.28972943 0.21795432 0.15881401 0.27302163
## [205] 0.25851175 0.41963736 0.25257437 0.32046549 0.29808320 0.21358326
## [211] 0.21854487 0.23244188 0.39786100 0.30210342 0.36501726 0.40908705
## [217] 0.30248886 0.34668628 0.23828245 0.46216300 0.34393431 0.46053080
## [223] 0.28432092 0.23436489 0.29730018 0.36149885 0.18062165 0.30389372
## [229] 0.35219163 0.34963635 0.43622642 0.30355318 0.35273627 0.22201769
## [235] 0.38025290 0.38068854 0.44682074 0.36889698 0.39142484 0.40025816
## [241] 0.16190456 0.38507810 0.43297728 0.44722257 0.36931782 0.37499346
## [247] 0.36716904 0.48399100 0.46192338 0.43163755 0.34416300 0.31644159
## [253] 0.33877788 0.47669515 0.34040470 0.38529433 0.15899874 0.42879516
## [259] 0.39489537 0.35930874 0.38273391 0.29719283 0.32917309 0.45405211
## [265] 0.39897218 0.40574637 0.42439340 0.46580989 0.39528830 0.45822270
## [271] 0.40751802 0.44533894 0.41713480 0.35047726 0.44441515 0.17916066
## [277] 0.17573744 0.28136123 0.41414124 0.35221648 0.35250881 0.44096824
## [283] 0.47369789 0.34587354 0.46538129 0.43125889 0.43772893 0.43172627
## [289] 0.40060880 0.39947123 0.41417213 0.18465626 0.41271631 0.40771847
## [295] 0.33290767 0.30341250 0.22719305 0.24123397 0.24147027 0.18777987
## [301] 0.34324898 0.18145060 0.01998837 0.06842105 0.17896670 0.16981424
## [307] 0.20749519 0.29168717 0.33090175 0.26634863 0.42319409 0.36871860
## [313] 0.39260848 0.36731836 0.35333230 0.28119576 0.40598159 0.42869998
## [319] 0.30438575 0.30269752 0.30717676 0.40593458 0.38441224 0.35750472
## [325] 0.41535187 0.30592744 0.41153293 0.22883304 0.30766150 0.41660286
## [331] 0.37407911 0.30118891 0.32037896 0.33130806 0.41528648 0.38319810
## [337] 0.34651253 0.28772042 0.35580009 0.43923832 0.47278158 0.25759640
## [343] 0.43932813 0.35151196 0.12203891 0.27728293 0.09017963 0.36911877
## [349] 0.34072943 0.41359510 0.43583626 0.30146271 0.39714938 0.37014761
## [355] 0.36893702 0.37350458 0.33459528 0.44473371 0.27131667 0.39849508
## [361] 0.43445562 0.33122323 0.38065764 0.33870815 0.23577121 0.23116024
## [367] 0.35514684 0.25881555 0.13985015 0.39492692 0.39304662 0.38227829
## [373] 0.41291584 0.34346882 0.44743770 0.44412297 0.37784581 0.45757765
## [379] 0.43836155 0.34870919 0.46780812 0.42088709 0.36905580 0.40609177
## [385] 0.27819558 0.32212463 0.18915709 0.38580468 0.38366454 0.33901819
## [391] 0.41109962 0.39241484 0.43072545 0.32779764 0.35903306 0.31233978
## [397] 0.42670875 0.44200466 0.38234083 0.44512775 0.37226788 0.25704286
## [403] 0.39510942 0.47167183 0.42541474 0.28667611 0.27362486 0.43007865
## [409] 0.45044834 0.36215556 0.34808447 0.39325444 0.38264619 0.28270202
## [415] 0.40820251 0.22353052 0.40917864 0.31635655 0.42887007 0.44959638
## [421] 0.34432747 0.38725935 0.39151862 0.35458826 0.38551508 0.38392397
## [427] 0.39052750 0.39760363 0.45568307 0.41898241 0.32150688 0.21843251
## [433] 0.29412492 0.24780440 0.26093338 0.30064206 0.43697931 0.28369380
## [439] 0.40279873 0.34819846 0.35895810 0.44549631 0.46108267 0.43944254
## [445] 0.35745147 0.41792324 0.27825146 0.27906687 0.23625013 0.23816053
## [451] 0.39051358 0.39460132 0.41489243 0.36201391 0.32071315 0.15815171
## [457] 0.32114754 0.42794379 0.26403013 0.36714083 0.34443393 0.42420938
## [463] 0.40323193 0.41190779 0.38134509 0.33840194 0.24470064 0.32303330
## [469] 0.26414014 0.44273993 0.23775858 0.31759709 0.29478737 0.42495425
## [475] 0.31653011 0.32111603 0.32568942 0.37932620 0.43512079 0.31921569
## [481] 0.32505583 0.29433855 0.25581681 0.21347212 0.28611138 0.37245965
## [487] 0.15863010 0.08161957 0.16639359 0.37313710 0.29871496 0.21645910
## [493] 0.32443687 0.11000024 0.10096446 0.34127065 0.44801744 0.22177817
## [499] 0.26014628 0.31129004 0.33880397 0.32574475 0.42678263 0.40940831
## [505] 0.46561616 0.37891909 0.28836471 0.31480402 0.35394206 0.42744294
## [511] 0.41121035 0.42295145 0.41996822 0.40697564 0.23355871 0.15322171
## [517] 0.32804203 0.25811799 0.20124368 0.34204520 0.37532332 0.35746346
## [523] 0.37117115 0.45226019 0.32721615 0.43273065 0.35227304 0.35386807
## [529] 0.37846667 0.41886099 0.40789170 0.35151333 0.46025638 0.36761879
## [535] 0.30891802 0.40750915 0.36779665 0.12166449 0.17264630 0.44656139
## [541] 0.41535553 0.25214437 0.21608011 0.37933563 0.40422222 0.44157381
## [547] 0.39726625 0.41678970 0.38361768 0.31224234 0.15133670 0.28879479
## [553] 0.40490784 0.42045345 0.41538779 0.37474210 0.33551199 0.36666548
## [559] 0.28761264 0.41113153 0.46649180 0.42780563 0.25296799 0.32228081
## [565] 0.25885574 0.36231627 0.42491581 0.35375603 0.42169115 0.22724126
## [571] 0.34989408 0.38950001 0.26741390 0.46689813 0.32763325 0.18357478
## [577] 0.27604038 0.23645678 0.43435755 0.36167834 0.24751073 0.28809842
## [583] 0.35705675 0.15614023 0.14861972 0.27578243 0.42067261 0.31096534
## [589] 0.31007638 0.19032399 0.34760922 0.38656124 0.37479418 0.36929867
## [595] 0.32910143 0.44182009 0.41855214 0.45564994 0.27729051 0.29591494
## [601] 0.45211424 0.21199559 0.45091145 0.47175252 0.37152178 0.17830592
## [607] 0.34957269 0.29008047 0.18645333 0.43759946 0.32046827 0.41558789
## [613] 0.36946811 0.33040838 0.48396684 0.35546787 0.26141879 0.40239317
## [619] 0.30877470 0.34848940 0.43964307 0.35096750 0.26150028 0.36365327
## [625] 0.41547087 0.38936249 0.41675380 0.40775087 0.31870839 0.30277155
## [631] 0.33272629 0.45230645 0.34314977 0.43409694 0.38108865 0.39734517
## [637] 0.32871616 0.44202428 0.34909893 0.37663075 0.32841963 0.27813086
## [643] 0.37256396 0.40768173 0.47285327 0.44252696 0.41744194 0.38079513
## [649] 0.41283334 0.46535355 0.38421374 0.30228087 0.44441046 0.39579502
## [655] 0.46804696 0.40919720 0.34898568 0.36075236 0.41610179 0.47119037
## [661] 0.36057145 0.31954097 0.36017224 0.43455562 0.26440106 0.11365672
## [667] 0.13927112 0.21512803 0.46413653 0.47222084 0.31020938 0.41529836
## [673] 0.35235405 0.45852756 0.46523903 0.25114708 0.43219760 0.46230308
## [679] 0.20279283 0.32297808 0.27238768 0.19494298 0.39049020 0.42831713
## [685] 0.47344185 0.37372088 0.25316878 0.23619191 0.39895604 0.28957250
## [691] 0.45377891 0.38543738 0.24058803 0.37609762 0.33029390 0.41368589
## [697] 0.39216522 0.23903535 0.26549492 0.45763102 0.39366331 0.11584876
## [703] 0.06202685 0.31203279 0.29228729 0.27377753 0.17709141 0.21454316
## [709] 0.12558162 0.12554666 0.21945617 0.14119521 0.13677408 0.04968259
## [715] 0.13785414 0.21296764 0.40412352 0.31697700 0.32771095 0.28736229
## [721] 0.46610302 0.39230093 0.44088115 0.29959472 0.02346418 0.12744407
## [727] 0.01867653 0.14634633 0.19653525 0.31087951 0.31027606
##
## $outpoint
## [1] 1 86 88 154 192 252 253 280 284 320 344 402 414 450 461 538 539 581 608
## [20] 609
##
## $medcurve
## [1] 248## $depth
## [1] 0.02020056 0.22785799 0.29899538 0.46275524 0.32301463 0.34794294
## [7] 0.07981302 0.10777832 0.18078641 0.17624304 0.33429730 0.22159675
## [13] 0.13287614 0.39191307 0.43137517 0.20434877 0.15878573 0.27419358
## [19] 0.37410257 0.18904547 0.30748012 0.41784497 0.46027051 0.42551609
## [25] 0.43695028 0.11572846 0.18655211 0.43883633 0.35169650 0.26609519
## [31] 0.40955511 0.43097506 0.47804832 0.46212465 0.40513609 0.39474991
## [37] 0.35940804 0.26778345 0.19758909 0.25930597 0.23305368 0.41792268
## [43] 0.41780317 0.24753177 0.34151341 0.39539125 0.43334272 0.35156202
## [49] 0.29749285 0.45909935 0.28266138 0.21485852 0.20341324 0.37436770
## [55] 0.46369589 0.45181077 0.26872852 0.13368969 0.15529030 0.26897377
## [61] 0.34325738 0.43275878 0.42862487 0.37890565 0.37962856 0.41522923
## [67] 0.31372162 0.36639883 0.42584403 0.40943539 0.43385529 0.46546019
## [73] 0.42999271 0.35141911 0.43734816 0.45991810 0.30643723 0.44645983
## [79] 0.39414526 0.40064671 0.31449413 0.26905539 0.33345491 0.37600087
## [85] 0.22901970 0.13512028 0.21075807 0.04383317 0.08494683 0.32635150
## [91] 0.21936667 0.33405746 0.21053650 0.44232832 0.18031083 0.31325371
## [97] 0.34191848 0.42101297 0.31621423 0.29865844 0.33843158 0.42449383
## [103] 0.41209011 0.37248183 0.25821674 0.29591479 0.40303629 0.36287096
## [109] 0.40380517 0.37864301 0.38960896 0.37688278 0.43878339 0.37716035
## [115] 0.24273043 0.30113485 0.40270121 0.40273145 0.37424513 0.32525911
## [121] 0.29993682 0.36776066 0.11752675 0.07402901 0.41224411 0.36964311
## [127] 0.20438554 0.40123440 0.39053392 0.29025235 0.30743644 0.33979586
## [133] 0.40248220 0.39556585 0.45387694 0.22975886 0.25125893 0.37877136
## [139] 0.45906091 0.43980906 0.45762286 0.40864507 0.43281147 0.46834058
## [145] 0.33937472 0.30107337 0.31145005 0.43326950 0.28750692 0.24668314
## [151] 0.08449957 0.30537813 0.35618304 0.26261938 0.28035655 0.36638172
## [157] 0.39066168 0.33911838 0.33234641 0.29003706 0.41455635 0.37976493
## [163] 0.44925713 0.43480054 0.26634756 0.30511188 0.42399177 0.45351878
## [169] 0.32321208 0.28023010 0.34656259 0.36456554 0.45902503 0.43271741
## [175] 0.32498629 0.45924170 0.44405864 0.42122927 0.32972266 0.42669614
## [181] 0.38441691 0.31650899 0.35755475 0.38247626 0.45784458 0.25273060
## [187] 0.28145525 0.34774872 0.44016047 0.33039784 0.33107736 0.33390693
## [193] 0.43291225 0.39277012 0.28776263 0.30068229 0.43938747 0.34253622
## [199] 0.22385466 0.33658962 0.43803662 0.44939775 0.39295600 0.32770787
## [205] 0.39696656 0.12014900 0.17814022 0.40893800 0.40425514 0.28174322
## [211] 0.45876271 0.35417083 0.39033398 0.29027873 0.38207268 0.48324057
## [217] 0.21687381 0.33976314 0.40703795 0.44220868 0.22868031 0.39527336
## [223] 0.38612117 0.40670382 0.39772520 0.21034687 0.33839043 0.47801570
## [229] 0.47052853 0.46591959 0.41554211 0.35722990 0.36895938 0.37293534
## [235] 0.37796180 0.31092026 0.42371034 0.44712428 0.23886577 0.38520976
## [241] 0.38411091 0.30743570 0.37721738 0.43029677 0.45564525 0.36883352
## [247] 0.35882674 0.32827855 0.28973751 0.34761131 0.28743964 0.36543042
## [253] 0.36102235 0.36538926 0.33585457 0.38400111 0.37184655 0.45269894
## [259] 0.32688942 0.19643881 0.40481737 0.45381160 0.38975320 0.42240812
## [265] 0.39079660 0.30466988 0.24885946 0.31067423 0.30460538 0.35100254
## [271] 0.26565513 0.38419356 0.32409117 0.17252355 0.31017805 0.42536788
## [277] 0.40702202 0.42483595 0.22705316 0.18603493 0.36287032 0.41631553
## [283] 0.36740200 0.36150945 0.44552637 0.44678321 0.46077764 0.44513176
## [289] 0.31793680 0.21727486 0.26900279 0.47217015 0.30325584 0.26120282
## [295] 0.38760494 0.25279525 0.46085149 0.39900324 0.39322295 0.46933082
## [301] 0.41352949 0.28206086 0.12913787 0.32654517 0.29016397 0.27624722
## [307] 0.44108416 0.44318845 0.28494546 0.34213772 0.16727740 0.30048120
## [313] 0.37430266 0.38857327 0.15047709 0.24984024 0.09723962 0.33337038
## [319] 0.16345953 0.10550021 0.32644490 0.29586132 0.21965073 0.37823140
## [325] 0.36275270 0.37787342 0.28193236 0.25032994 0.33504180 0.17525760
## [331] 0.06060683 0.20822808 0.26465977 0.34597703 0.33780524 0.23422155
## [337] 0.02869530 0.28414989 0.39215824 0.36008723 0.37604963 0.30152794
## [343] 0.34765706 0.14410307 0.32555571 0.35376794 0.38990419 0.35222388
## [349] 0.39903727 0.44996624 0.30619243 0.35407727 0.30302800 0.09006114
## [355] 0.28491551 0.23067672 0.34277275 0.16204119 0.29120377 0.42125075
## [361] 0.42773346 0.22649101 0.28662501 0.36335933 0.30515793 0.30071320
## [367] 0.43408803 0.24053837 0.20448306 0.25717697 0.42219960 0.38655385
## [373] 0.37970199 0.41888010 0.14259712 0.36430152 0.19369039 0.34136743
## [379] 0.30976808 0.30063828 0.36054581 0.05198785 0.15095986 0.37221521
## [385] 0.44337372 0.37215235 0.39918390 0.35844210 0.16235594 0.11468261
## [391] 0.30918731 0.05948398 0.20741586 0.25519892 0.04068611 0.13061132
## [397] 0.32409803 0.27047246 0.30017356 0.11739439 0.17673662 0.04481809
## [403] 0.30991437 0.34857221 0.31116020 0.42128741 0.32183081 0.43416941
## [409] 0.37748274 0.30055942 0.42404885 0.41673421 0.33606346 0.14464051
## [415] 0.08964509 0.36755934 0.32578568 0.10324950 0.28520829 0.39071931
## [421] 0.39031431 0.39382685 0.28871823 0.20562332 0.09911566 0.32338325
## [427] 0.29265073 0.41477981 0.40058174 0.42696458 0.43951515 0.35245161
## [433] 0.23882046 0.20815467 0.16056606 0.38719694 0.30625143 0.39933603
## [439] 0.30891694 0.37217551 0.42359508 0.47689886 0.42067490 0.32746169
## [445] 0.17494625 0.17480681 0.45715626 0.39988816 0.23225853 0.12484777
## [451] 0.15310768 0.28611188 0.41813052 0.31900290 0.36147907 0.46572915
## [457] 0.38166720 0.46044163 0.43309687 0.19028982 0.33241873 0.36762761
## [463] 0.41920588 0.29948537 0.30874639 0.41127399 0.21651357 0.32065548
## [469] 0.30479976 0.40226125 0.31797890 0.46315499 0.36953779 0.34838567
## [475] 0.32816879 0.39474987 0.38889630 0.41955205 0.39320339 0.32748592
## [481] 0.36646696 0.35847429 0.25752215 0.40209983 0.37733138 0.34534437
## [487] 0.44746850 0.36849016 0.40159527 0.39986942 0.41124883 0.44846306
## [493] 0.42012185 0.45413403 0.33332009 0.43285699 0.42000025 0.42590091
## [499] 0.37424158 0.27749787 0.40735181 0.38979561 0.36262226 0.38126477
## [505] 0.33628822 0.30055367 0.45259965 0.41546326 0.32108037 0.26423044
## [511] 0.46940006 0.40018895 0.47639120 0.35841738 0.30034731 0.43772209
## [517] 0.39678922 0.38626143 0.28511769 0.42757512 0.47496973 0.24871418
## [523] 0.20261685 0.38829756 0.33760628 0.21572373 0.39440262 0.43833114
## [529] 0.45459814 0.46129546 0.34632603 0.38141591 0.37052205 0.44125130
## [535] 0.37711842 0.47545195 0.36319740 0.33039116 0.17968726 0.39094589
## [541] 0.39183248 0.45175699 0.32918187 0.29415639 0.37537805 0.31801928
## [547] 0.24890838 0.45669301 0.42579453 0.37240941 0.31707618 0.47057079
## [553] 0.27349306 0.35546441 0.29249691 0.33433473 0.47615656 0.47603106
## [559] 0.44174836 0.41500098 0.37667749 0.44815845 0.23532387 0.44742537
## [565] 0.37771199 0.31375439 0.31031473 0.30360627 0.45580150 0.42158824
## [571] 0.43502879 0.39129782 0.37465085 0.43022434 0.32652367 0.36139759
## [577] 0.45652895 0.43502249 0.22318546 0.27299696 0.41991527 0.41822522
## [583] 0.21387785 0.36928929 0.45115818 0.44372708 0.35273169 0.35275786
## [589] 0.36710136 0.30012856 0.35337278 0.34832857 0.38840068 0.38700516
## [595] 0.39525399 0.18124554 0.20886200 0.34781369 0.34381408 0.20054564
## [601] 0.25122638 0.31626351 0.32375635 0.41331508 0.25773484 0.29254938
## [607] 0.43096205 0.37749867 0.22397270 0.29618902 0.33308839 0.32890824
## [613] 0.34295050 0.33617193 0.43918056 0.47977660 0.29946035 0.39900093
## [619] 0.44030845 0.34173877 0.33234853 0.41928596 0.36905528 0.26070751
## [625] 0.22134241 0.18177196 0.30305724 0.42458769 0.37474217 0.42290918
## [631] 0.37038876 0.35784857 0.44646446 0.31705754 0.23455875 0.40558413
## [637] 0.29706243 0.34985942 0.17652669 0.16039927 0.18181460 0.36741819
## [643] 0.36072242 0.30645022 0.35512851 0.28093834 0.23633225 0.40059358
## [649] 0.39662009 0.33081860 0.27872771 0.42799492 0.45703214 0.23502713
## [655] 0.39404879 0.46912870 0.39224482 0.23906363 0.36315045 0.36313552
## [661] 0.10794215 0.40656805 0.41933838 0.27047958 0.07678777 0.11941808
## [667] 0.18016191 0.28373036 0.33236204 0.34154147 0.32979702 0.45686912
## [673] 0.39084321 0.34518121 0.39053559 0.36412840 0.43903197 0.29443433
## [679] 0.31246238 0.26874563 0.43035615 0.39452581 0.29195150 0.29946614
## [685] 0.34819196 0.41356890 0.25672926 0.44048331 0.28527606 0.26322428
## [691] 0.31178193 0.44363419 0.22390119 0.35965507 0.36311026 0.18427059
## [697] 0.38824756 0.41045538 0.21017307 0.31493605 0.32359840 0.20175499
## [703] 0.19705314 0.43298054 0.45634597 0.44562222 0.40562009 0.46482185
## [709] 0.34758967 0.33397561 0.35044933 0.41891298 0.24813091 0.18794054
## [715] 0.19099057 0.22557162 0.40735258 0.37654295 0.45078187 0.30098379
## [721] 0.21210699 0.15988473 0.31262529 0.23704439 0.23987478 0.16277530
## [727] 0.32040125 0.44751366 0.46662765 0.39010680 0.37643819
##
## $outpoint
## integer(0)
##
## $medcurve
## [1] 216
par(mfrow=c(1,1))
# Depth calculation example for 'cnt'
library(fda.usc)
cnt_fdata <- fdata(smoothed_fd_list$cnt)
cnt_fd_depth <- depth.FM(cnt_fdata, trim = 0.25)
depth_values <- cnt_fd_depth$dep
summary(depth_values)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.1048 0.3891 0.4868 0.5000 0.6182 0.8584Depth Comparisons and Wilcoxon Tests
# Compute holiday and season labels as before
holiday_labels <- bike_data %>%
group_by(day) %>%
summarise(is_holiday = first(is_holiday)) %>%
arrange(day) %>%
pull(is_holiday)
wilcox_test_holiday <- wilcox.test(depth_values[holiday_labels==1], depth_values[holiday_labels==0])
season_labels <- bike_data %>%
group_by(day) %>%
summarise(season = first(season_label)) %>%
arrange(day) %>%
pull(season)
wilcox_test_season <- wilcox.test(depth_values[season_labels=="Summer"], depth_values[season_labels=="Winter"])
cat("Wilcoxon Test for Holidays vs Non-Holidays (cnt):\n")## Wilcoxon Test for Holidays vs Non-Holidays (cnt):##
## Wilcoxon rank sum test with continuity correction
##
## data: depth_values[holiday_labels == 1] and depth_values[holiday_labels == 0]
## W = 2757, p-value = 0.0003906
## alternative hypothesis: true location shift is not equal to 0##
## Wilcoxon Test for Summer vs Winter (cnt):##
## Wilcoxon rank sum test with continuity correction
##
## data: depth_values[season_labels == "Summer"] and depth_values[season_labels == "Winter"]
## W = 11814, p-value = 1.108e-06
## alternative hypothesis: true location shift is not equal to 0depth_df <- data.frame(
day = unique_days,
depth = depth_values,
holiday = factor(holiday_labels, levels = c(0,1), labels = c("Non-Holiday", "Holiday")),
season = factor(season_labels, levels = c("Spring", "Summer", "Fall", "Winter"))
)
p_depth_holiday <- ggplot(depth_df, aes(x = holiday, y = depth, fill = holiday)) +
geom_boxplot(alpha = 0.7) + theme_minimal() +
labs(title = "Depth by Holiday vs Non-Holiday", x = "Holiday Status", y = "Depth")
p_depth_season <- ggplot(depth_df, aes(x = season, y = depth, fill = season)) +
geom_boxplot(alpha = 0.7) + theme_minimal() +
labs(title = "Depth by Season", x = "Season", y = "Depth")
# include t1 and cnt in daily averages
daily_avgs <- bike_data %>%
group_by(day) %>%
summarise(
mean_cnt = mean(cnt, na.rm = TRUE),
mean_t1 = mean(t1, na.rm = TRUE),
mean_t2 = mean(t2, na.rm = TRUE),
mean_hum = mean(hum, na.rm = TRUE),
mean_wind = mean(wind_speed, na.rm = TRUE)
) %>%
arrange(day)
depth_df <- depth_df %>%
mutate(
mean_t1 = daily_avgs$mean_t1,
mean_t2 = daily_avgs$mean_t2,
mean_hum = daily_avgs$mean_hum,
mean_wind = daily_avgs$mean_wind,
mean_cnt = daily_avgs$mean_cnt
)
p_depth_t2 <- ggplot(depth_df, aes(x = mean_t2, y = depth)) +
geom_point(alpha = 0.7, color = "darkgreen") + geom_smooth(method = "loess", se = TRUE, color="red") + theme_minimal() +
labs(title = "Depth vs Mean T2", x = "Mean T2 (°C)", y = "Depth")
p_depth_hum <- ggplot(depth_df, aes(x = mean_hum, y = depth)) +
geom_point(alpha = 0.7, color = "purple") + geom_smooth(method = "loess", se = TRUE, color="red") + theme_minimal() +
labs(title = "Depth vs Mean Humidity", x = "Mean Humidity (%)", y = "Depth")
p_depth_wind <- ggplot(depth_df, aes(x = mean_wind, y = depth)) +
geom_point(alpha = 0.7, color = "brown") + geom_smooth(method = "loess", se = TRUE, color="red") + theme_minimal() +
labs(title = "Depth vs Mean Wind Speed", x = "Mean Wind Speed (km/h)", y = "Depth")
# New plots for mean_t1 and mean_cnt
p_depth_t1 <- ggplot(depth_df, aes(x = mean_t1, y = depth)) +
geom_point(alpha = 0.7, color = "blue") + geom_smooth(method = "loess", se = TRUE, color="red") + theme_minimal() +
labs(title = "Depth vs Mean T1", x = "Mean T1 (°C)", y = "Depth")
p_depth_cnt <- ggplot(depth_df, aes(x = mean_cnt, y = depth)) +
geom_point(alpha = 0.7, color = "darkorange") + geom_smooth(method = "loess", se = TRUE, color="red") + theme_minimal() +
labs(title = "Depth vs Mean Cnt", x = "Mean Cnt", y = "Depth")
# have 7 plots in total:
# 1. p_depth_holiday
# 2. p_depth_season
# 3. p_depth_t2
# 4. p_depth_hum
# 5. p_depth_wind
# 6. p_depth_t1
# 7. p_depth_cnt
# Arrange them in a suitable grid. For 7 plots, do a 3x3 grid with some dummy plots:
dummy_plot1 <- ggplot() + theme_void() + labs(title="")
dummy_plot2 <- ggplot() + theme_void() + labs(title="")
grid.arrange(
p_depth_holiday, p_depth_season, p_depth_t2,
p_depth_hum, p_depth_wind, p_depth_t1,
p_depth_cnt, dummy_plot1, dummy_plot2,
ncol=3
)
Deepest Curves for All Variables
all_variables_for_depth <- c("cnt", "t1", "t2", "hum", "wind_speed")
par(mfrow=c(2,3))
for (var in all_variables_for_depth) {
if (var %in% names(results) || var=="t1") {
var_data <- hourly_data %>%
select(matches(paste0("^", var))) %>%
as.matrix()
var_data[is.na(var_data)] <- 0
var_fd <- smooth.basis(seq(0,23,length.out=ncol(var_data)), t(var_data), basis_depth)$fd
var_fdata <- fdata(var_fd)
var_depth <- depth.FM(var_fdata, trim=0.25)$dep
ordered_depth_indices <- order(var_depth, decreasing = TRUE)
cat("\nTop 5 Deepest Curves for:", var, "\n")
print(ordered_depth_indices[1:5])
plot(var_fdata[ordered_depth_indices[1:5]],
main = paste("Top 5 Deepest -", var),
col = rainbow(5))
} else {
cat("\nVariable", var, "not found. Skipping.\n")
}
}##
## Top 5 Deepest Curves for: cnt
## [1] 454 455 678 69 517##
## Top 5 Deepest Curves for: t1
## [1] 263 655 124 465 331##
## Top 5 Deepest Curves for: t2
## [1] 263 704 655 331 650##
## Top 5 Deepest Curves for: hum
## [1] 615 248 670 254 685##
## Top 5 Deepest Curves for: wind_speed
## [1] 216 513 616 33 557
Clustering Based on FPCA Scores
# Assume that smoothed_fd_list has fd objects for: cnt, t1, t2, hum, wind_speed
# Compute depth values for each variable and store in depth_values_by_var
depth_values_by_var <- list()
depth_vars <- c("cnt", "t1", "t2", "hum", "wind_speed")
for (v in depth_vars) {
var_fdata <- fdata(smoothed_fd_list[[v]])
var_fd_depth <- depth.FM(var_fdata, trim = 0.25)
depth_values_by_var[[v]] <- var_fd_depth$dep
}
# Clustering Based on FPCA Scores
vars <- c("cnt", "t1", "t2", "hum", "wind_speed")
# depth_values_by_var is defined and contains depth values for all variables
cluster_plots <- list()
set.seed(123)
k <- 3 # number of clusters
for (v in vars) {
# Extract FPCA object and dynamic_nharm for the variable
fPCA_obj <- pca_results[[v]]$fPCA
dynamic_nharm <- pca_results[[v]]$dynamic_nharm
# Extract scores for the chosen PCs
scores_var <- fPCA_obj$scores[, 1:dynamic_nharm]
# Run k-means clustering
km_res_var <- kmeans(scores_var, centers = k, nstart = 10)
# Create a depth dataframe for this variable
depth_df_var <- data.frame(
day = unique_days,
depth = depth_values_by_var[[v]],
cluster = factor(km_res_var$cluster)
)
# Create boxplot for this variable
p_cluster_var <- ggplot(depth_df_var, aes(x = cluster, y = depth, fill = cluster)) +
geom_boxplot(alpha = 0.7) + theme_minimal() +
labs(title = paste("Depth Distribution by Cluster -", v), x = "Cluster", y = "Depth")
cluster_plots[[v]] <- p_cluster_var
}
# Arrange all 5 cluster plots in a 2x3 grid with a dummy plot if needed
if (length(vars) < 6) {
dummy_plot <- ggplot() + theme_void() + labs(title = "")
cluster_plots <- c(cluster_plots, list(dummy_plot))
}
grid.arrange(grobs = cluster_plots, nrow = 2, ncol = 3)
Conclusion
This comprehensive FDA analysis provided insights into patterns, effects of season/holiday/weather, and underlying functional structure. Depth analysis and clustering further enriched the understanding.
Appendix
Metadata Description
cat("
**Metadata:**
- **timestamp**: Timestamp.
- **cnt**: The count of new bike shares.
- **t1**: Real temperature in °C.
- **t2**: Temperature in °C 'feels like'.
- **hum**: Humidity in percentage.
- **wind_speed**: Wind speed in km/h.
- **weather_code**: Weather category.
- **is_holiday**: 1 if holiday, 0 otherwise.
- **is_weekend**: 1 if weekend, 0 otherwise.
- **season**:
- 0: Spring
- 1: Summer
- 2: Fall
- 3: Winter
**weather_code** categories:
- 1 = Clear
- 2 = Scattered Clouds
- 3 = Broken Clouds
- 4 = Cloudy
- 7 = Rain
- 10 = Rain with Thunderstorm
- 26 = Snowfall
- 94 = Freezing Fog
")Metadata:
- timestamp: Timestamp.
- cnt: The count of new bike shares.
- t1: Real temperature in °C.
- t2: Temperature in °C ‘feels like’.
- hum: Humidity in percentage.
- wind_speed: Wind speed in km/h.
- weather_code: Weather category.
- is_holiday: 1 if holiday, 0 otherwise.
- is_weekend: 1 if weekend, 0 otherwise.
- season:
- 0: Spring
- 1: Summer
- 2: Fall
- 3: Winter
weather_code categories: - 1 = Clear - 2 = Scattered Clouds - 3 = Broken Clouds - 4 = Cloudy - 7 = Rain - 10 = Rain with Thunderstorm - 26 = Snowfall - 94 = Freezing Fog